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ABSTRACT 

Mathematical models are used to represent a 
physical system hut for most of the real systems, it is 
very rare, that a mathematical model represents all the 
details of the actual process. The models that are based 
on physical and chemical principles such as conservation 
of mass, momentum and energy, chemical kinetic mechanisms 
etc. are, sometimes, called as mechanistic models. Parameter 
estimation forms an integral part of model building for such 
cases. 

Parameter estimation procedures for algebraic 
models are quite well established but insignificant literature 
exists for the case of differential equation models. In the 
present work, four different algorithms for parameter estima- 
tion in differential models are compared for their computational 
®ase and the reliability of the parameter estimates. Approximate 
analytical solutions are also presented for models linear in 
state variables. 

Two models representing batch reactor kinetic 
data with eight data sots representing different experimental 
conditions, were used to test the stability of the algorithms. 

Two of the algorithms employed explicit and 
implicit method of integration respectively, combined with 
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nonlinear least squares analysis. The third algorithm 

( 29 ) 

used the Constrained Simplex method for minimization 

of an objective function based on variance- covariance matrix 
of residuals. A weighted residual method was used in the 
fourth algorithm which gave explicit estimates of parameters, 
thereby, obviatir^the necessity of numerical integration and 
function minimization. 

It is observed that the weighted residual method 
(Algorithm 4) is quite satisfactory for the task of parameter 
estimation in differential equation models as the parameters 
obtained from the application of this method had the m«e*imum. 
reliability and minimum computer central processing time. It 
is hoped that this work would be helpful to the future inves- 
tigators working with dynamic models. 



CHAPTER 1 


IETROKJCTIOE 

One of the important aspects of Chemical 
Engineering systems is. their overall design and analysis. 
These complex systems are represented by mathematical models 
which may be derived in a variety of ways, for example, 

a. They may be specified directly from basic principles 
of Physics and Chemistry i.e. knowing the underlying phy- 
sical and chemical principles of the system one can write 
down the model equations. Models of this type, relevant to 

chemical engineering systems, are often called " transport 

( 1 ) 

phenomena models" . Models given by the phenomenological 
equations of change, that is the continuum equations descri- 
bing the conservation of mass, momentum, and energy, are 
examples of transport phenomena models. These models may 
also be termed as deterministic models as the variables and 
parameters involved can be assigned a definite fixed number, 
for any given set of conditions. 

b. At the other extreme is the so called " Black Box" model 
where essentially no a-priori information is available about 
the system. The model must be determined solely on the basis 
of the experimental data gathered. Polynomial expressions 
used to fit the data collected from an experiment represent 

a class of empirical or " Black Box" models. There is a lot 
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of 'uncertainty about these models and extreme caution, there- 
fore, should be employed while using them for the prediction 
of the state of the system for some other conditions not used 
at the time of original experimentation. 

c. Most real situations lie somewhat intermediate of these 
two extremes. We often have some information on the basic 
nature of the system but not enough to set down precisely 
the relation among all the variables. It is in this case 
that the modeling procedure often includes extensive comp- 
arison of alternate models that can be written based on diff- 
erent underlying mechanisms, to determine the ! L best model ” . 
This is known as mechanistic model building. It is also in 
this case that the stochastic models find extensive usefulness. 
In a typical Chemical Engineering situation, for example, there 
might exist a mathematical model n =f (9,^ ) that is perhaps 
nonlinear in the parameters 9, which relates a measurable res- 
ponse n to the controllable variables S, . The response might 
be the rate of a chemical reaction, the variables might be 
partial pressures of the reactants and products, and the para- 
meters might be adsorption equilibrium constants and an overall 
reaction rate constant. we may not know the true values of the 
parameters involved in the model and this introduces some uncer- 
tainty in the mechanistic models. In some way this uncertainty 
is minimized to give a reasonably accurate model for the system 
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under consideration, which may he used with some degree of 
confidence for extrapolation. 

Once the mod. els are formulated the next impor- 
tant thing is the solution of these models to indicate the 
nature of the response of the systems concerned. This may 
require the use of analytical or numerical techniques. Once 
the response is predicted from a model it may he compared 
with the experimental data to find out the most plausible 
model representing the given situation. The techniques inv- 
olved in the solution of the models would depend upon their 
types, for example, the relationship between the controllable 
variables and the dependent variable or response might be rep_ 
resented by algebraic equations-meaning that the response of 
the system is measured experimentally. In such situations, . 
the comparison between the model prediction and the data can 
be carried out directly based on some standard statistical 
techniques such as linear and nonlinear least squares analysis. 
These models may be of two different types depending upon the 
complexity of the system studied. Tor example, if the system 
considered represents a single reaction between some chemical 
compounds, carried out in the presence of a catalyst, then the 
model representing this physical situation would be a " single 
response model", as the rate measurement of only one component 
is sufficient to identify the system completely. Take, for 
example, the catalytic dehydrogenation of ethyl alcohol into 
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acetaldehyde and hydrogen. Depending on whether 1) the 
adsorption of ethyl alcohol, A, 2) the surface dehydrogen- 
ation, or 3) the desorption of acetaldehyde, S, is the 

rate determining step, the following mechanistic models of 

( 2 ) 

the Hougen-Watson type have "been derived 7 , in terms of 
equilibrium constants, rate constants and partial pressures, 
c 2 h 5 oh ir— CH^CHO + h 2 (1) 

A R 

E(r A ) = k(p A -p2 / K )/ (1+ E A p R 2 / K+(K R +K s )p R ) (2) 

E(r A ) = kK A (p A -p2 /£)/ (1 +K a p a + (K R +K s )p R ) 2 (3) 

S(r A ) = K(p A -p 2 /K)/ (p R+ K A P APit+ KK a p A+ K sPR ) 2 (4) 

where E (r A ) = expected value of the rate of ethanol dehydro- 
genation. 

To illustrate the multiresponse situation, we 
may consider the following moderately complex structure of 
solid catalyzed chemical reactions . 

' D 

Among many others, three representative rival mechanisms may 
be thought of: 

1) all reactions are described by simple first order homogene - 
ous kinetics i.e. , interaction of the species with the catalyst 
is not explicity accounted for; 
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2) it is supposed that the reaction A^—^-'B+C occurs on dual 
sites, whereas A”~ D and C — *-E both only require a single 
site; this model is the true model. 

3) all reactions occur on single sites, whereby e.g. the 
component B is not adsorbed noticeably on the catalyst 
surface. 

The following rate equations are then derived 
( table 1). The different reaction rates are the responses 
of the system. 

Although there are five components involved, 
there are only three independent chemical reactions. 


Model I 

Model II 

r D = k 1 (P A -PD /K) /Pt 

r D = k 1 K A (P A -P D / K ^^ 

r B = k 2 (P A-PB P C /K,)/P t 2 

r B= k 2 K A {v A- v B I, C /K '' )/i ‘ 2 

r C = k 3 p c/ p t 

r c= k 3 K e> p c /0 

0 = 1 +K A P A +K bPb +K cPc +E D P D +K ePe 

Model III 

Model IV 

r D“ k 1 K A^ P A -P D /K) /^ 

r D = k 1 K A (Pa-PiAW 

r B = k 2 K A^ Pa“ p bPc /K ' ) 

r B= k 2 K A ( P A ~PbPc/ K ' ) / V-' / Pt 

r c = k 3 K c P c / H' 

r c = k 5 K c P c /’r' 


^ = 1 + e aPa +k cPc+ k dPd+ k ePe 


Table 1. Rival Rate Equations for A^B +C — >E 

^ D 
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In some situations the equations may "be of 
differential type and involve the controllable variables 
implicitly. An example representing this situation, is the 
acid-catalyzed hydrolysis of propylene oxide which, in dilute 

,,.r- 

( 3 ) 

solution is described by the single reaction v ' . 

^ '■'**»* IT i , 

CH^CH CH 2 + H 2 0 -Si. CH 3 CH0HCH 2 0H (6) 


If the reaction is carried out adiabatically in a batch reactor, 
the walls of which have negligible heat capacitance, the sys- 


tern equations are 

of the 

form: 




dt 

- -®i 

CD 

•a 

k i +y 2 

■) y? 

(7) 

dy 2 

d¥~ 

= je 1 

CD 

H 

© 0 y 9 

M+y 2 ; 

y^ 

(8) 


with initial conditions 


y 1 = y 1 


t = 0 


( 9 ) 


y 2 = y 2o 

where y^ is the molar concentration of propylene oxide, y 2 is 
a normalised temperature, t is time; ©^and © 2 are undeter- 
mined parameters related to the pre-exponential Arrhenius 


factor k Q and activation energy A by 

©1 = k Q exp ) , © 2 = 


( 10 ) 


J = is a thermicity factor, measured 

s p d 

independently, and n is the reaction order, assumed to be 
known. 
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In tlie previous example of multiresponse system, 
the continuity equations for the species D,B and E in a plug 
flow integral tubular reactor read: 

dx— dXg dXg 

dw/I' Ao = r D’ dw/P Ao = r B 1 dw/P Ao r E 5 X JT X B =X E = 0 

( 11 ') 

when w/F Aq = ® 

When the expressions for r-, r_ and r_, as given in table 1, 
are substituted in Eq. (11), we get a set of simultaneous 
differential equations. 

As we see that it is not uncommon to have 
many rival mechanisms for the same system, due to the uncer- 
tainty in the theoretical knowledge of physical situations, 

it is worthwhile to have a generalised strategy for model 
discrimination for example in the case of heterogeneous 

kinetic systems. Such strategy . would involve a comparison of 
several rival kinetic models to select the most " adequate" 
model, and the estimation of parameters involved in this model. 
The method might depend upon the way the experimental data 
are gathered which, in turn, may depend on the way the exper- 
iments are planned. There are two ways in which experiment- 
ation can be planned so as to allow the maximum discrimination 
between the rival models - sequential or non sequential experi- 
mental designs. The non sequential way would mean the a-priori 
selection of experimental conditions in terms of controllable 
variables whereas sequential experimental design strategy 
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might be used to reduce the experimental effort required at 
the same time carrying out the discrimination and parameter 
estimation sequentially. However, this procedure would requ- 
ire a knowledge of a few preliminary data points which may 
be obtained by using nonsequential fractional factorial desi- 
gns. 

A typical example for the case of model discri- 
mination may arise when an experimenter is studying a chemical 
reaction in which a reactant A is used to make a product B. 
Simultaneously, however, an undesired by product C is formed. 
Suppose further that the experimentor knows that the reaction 
is one of the following: 

Mechanism 1 : A — B — * 0 

( 12 ) 

Mechanism 2 : A —*■ B - r -r>- C 

If one of these mechanisms is actually correct 
and all experiments consists of measurements of the concentra- 
tions of B at different times, then, where to place the exper- 
iments to discriminate most efficiently between the two models 
will depend in general on the values of the estimated parame- 
ters. A sequential scheme therefore suggests itself in which 
at each stage the two models are refitted and a decision is 
made as to what experiment should be run next. 

Designs for parameter estimation are carried out 
if the form of the theoretical model is known, the problem 
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which confronts the experimentor is to evaluate the physical 
parameters (e.g., the rate constants in chemical kinetics). 

If the experiments are not carefully planned the experimental 
points may he so situated in the space of the variables that 
estimates that are obtained for the parameters are not only 
imprecise but also highly correlated. Once the data are 
collected a statistical analysis, no matter how elaborate , 
can do nothing to remedy this unfortunate situation. However, 
by selecting a suitable experimental design in advance ,these 
shortcomings can be overcome. The initial set of experiments 
could be run according to a standard fractional factorial 
pattern ^ ^ and then sequential strategy implemented based 
on the initial estimates obtained. 


A flow diagram that illustrates the sequential 


design procedure for model discrimination and parameter esti 
mation can be drawn as below: 


Re jectea 
Model(s) 


Initial set of 
experiments based 
on fractional 
factorial design 


Estimate parametersj 
for all the rival j 
models | 




\ 


''Some 


k.— ~ standard^ 
'•v^deq.uacy tes. 


Perform addit- 
ional experim- 
ent based on 
some discrimi- 
nation criter- 
ion. 


t yes 



Accepted 
models (s) 
If more 
than one? 


Ho 


Most 
adequate 
model 


Sequential Design Strategy for Model Discrimination 
and Parameter Estimation. 


Eig.1 
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A -considerable work has been done in the area of 

> ■ 

parameter estimation and model discrimination for the case of 
algebraic models v ' but meagre information seems to be avai- 
lable for the case of models governed by differential equations. 

The estimation of parameters in ordinary differe- 
ntial equations is of considerable importance , for example, in 
the determination of kinetic parameters from batch and integral 
reactor data, studying the dj^namics of systems such as missiles 
and spacecraft trajectories from tracking data, and in several 
other occasions where the rate of change of state variables 
w.r.t. some independent variable cannot be measured directly. 

The procedure, in many cases, is usually the extension of the 
principles used for the case of algebraic equations. However , 
the task of parameter estimation and model discrimination, in 
the case of the former which is quite complicated due to the 
multiresponse systems involved, is further complicated due to 
the nonavailability of the responses directly from model thus 
making it necessary either to integrate the differential equa- 
tions or to use some variational principles which may circum- 
vent the integration of the equations. A redefinition of 
dependent and independent variable would be appropriate at 
this stage. The dependent variable for the case of algebric 
equations is replaced by the derivative of the state variable 
in case of differential models; for example, it could be the 
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immeasurable rate of change of the concentration of reactants 
and products w.r.t. time or space -which, incidently, are the 
independent variables in differential equations. The concen- 
trations or such variables in which a change is observed w.r.-fc. 
independent variables, are known as state variables as they 
determine the state of the system at a given time or space. 

The experimental data that is collected, in this situation, 
may be a continuous function of time and / or space or might 
have been gathered only at discrete points. The statistical 
analysis of these two types of situations is based on concept 
of " 'noise" that gets added to the data during measurements 
due to the inaccuracy of measuring instruments and human 
oversight. This "noise" may be of two types as illustrated 
below: 

a. Noise in the output of the system; for example, random 
errors in the measurements of reaction rate or concentrations 
in a reaction kinetic study. 

b. Dynamical noise that gets added to the input of the 
system implying that the differential equations act on this 
noise too, along with the state variables. A pictorial 
representation would he: 


Error 


x(-b) 


determini- 
stic model 


(a) 


y( 






Y(t) X (*) 


De termini- ! Y(t ) 
stic model! — * 


(b) 



12 


where the bold letters denote random variables and stands 
for the error structure superimposed on the output of the 
system* 

We will deal only with such systems that are 
represented by ordinary differential equations and do not 
involve any dynamical noise, as these systems are quite comm- 
on in kinetic reaction studies of complex reaction networks. 

LITERATURE REVIEW : 

The advent of high speed digital computers has 
seen a large number of research papers published in the area of 
parameter estimation and model discriminination for the case of 
algebraic models but the literature on parameter estimation and 
model discrimination in multiresponse systems and differential 
equation models appears to be meagre. 

The problem of parameter estimation can be posed 

as follows: 

The system is described by a vector of ordinary 
differential equations 

k = £ (£,k ) 

( 13 ) 

x(o)= 4 

where I i« an n-dimensional state vector and k is the p- 
dimensional parameter vector. The output of the system is 
the m-dimensional vector Measurements of the output are 
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made at R times t^,t^# •»» t^. The measurements are related 
to the state vector according to 

Y r = h (t,I (t r ) + ( Measurement Errors) (14) 

r = 1,2,..., R„ 

Given the observations Y r we want to fit the model given 
by Eq.(l3) by the data Y r r = 1,2,.,,R in some optimal way* 

If we assume the initial state of process X Q to 

be known, we desire to select the parameters k to minimize 

the least squares criterion ( if X is also unknown we will 

want to select both X Q and k ) 

R 

s= jSjtlSr-fe ( t r , 2 (t r )) jj 2 Sr (15) 

where Q r is an mm positive semidefinite weighting matrix 
incorporating the relative accuracy of each component of the 
measurement 'vector at each time of measurement. 

There are essentially eight different ways by 
which we can estimate the parameters involved in heterogeneous 
reacting systems for differential types of models, from exper- 
imental data. These are discussed in the literature as follows: 

1 . Analytical ( exact or approximate ) integration of the 
set of differential equations, and subsequent application of 
iterative nonlinear least squares regression techniques (5-9) ^ 
Kittrell et.al have proposed a method, which is a varia- 
tion of the method of integration, that transforms the dependent 
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variable to achieve an error distribution as consistent as 
possible with the assumptions inherent in the least squares 
analysis. The method can be applied to the determination 
of unknown reaction order and rate constants for power latf 
rate models very successfully and for a small set of equations. 

2a* Differentiation of the emprical data directly, and subse- 
quent application of linear least squares regression techniques 
( 10*7» 1 1 , 1 2) ^ Here one has to bear in mind that as differentia- 
tion is more error prone than the original data and so the esti- 
mates obtained, may not be the optimal estimates of the rate 
constants. Thus, differentiation of empirical data increases 
the existing inherent errors in the experimental data and there- 
fore, category 2 is not satisfactory. 

2b. Another procedure that can be followed based on differen- 
tiation is to use linear regression to fit the emprical data, 

differentiation of the regression equation followed by linear 

fix') 

regression to estimate the coefficient v * 

3a. Numerical integration of the set of differential equations 
using emprical data directly, followed by iterative non-linear 

least squares technique (’•4-17), Ball and Groenweghe^ ^ bave 
used the values of the constants, obtained by the differentia- 
tion of the data and subsequent use of linear least squares 
equation, as the initial guess values for the iterative numeri- 
cal integration procedure. 
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5b. linear regression to fit the emprical data, followed by 
numerical integration of the set of differential equations. 


4. Trial and error search using analog computers to match 
the emprical data (^,19). 

5. Method of differential corrections may also be used . 
Incase, if the differential equations admit explicit solutions 
in terms of known elementary or transcendental functions then 
according to the method, if 


y = f ( t,-<, p ,i ) 


is the solution of the differential equation and if i°>, p {0) 
and are given estimates of the unknown parameter values, 

the n equations 

f (t,,q ,?,$)= f (t, 

ct f 


(o) *(o) . ( o) \ 

9 p t irX J 


m 


+ 2 . ~ (-< k - * k (o) ) 
k=1 3*k K * 


m , ' ( o) m % „ ( o) 

+ f=i sr < > + f =1 ^ > 


( 16 ) 


are formed where the derivatives are evaluated at (t * ^ 

The equations give the corrections to the first approximation 
and then, the new values of the estimates are used and the 
process continued iteratively. The method may be employed 
in the case of differential being nonlinear for which explicit 
solutions are not available . 
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( 2122 ') 

6. Quasilinearisation v * 1 is another technique which 

has been used to tac£Le the problem. In a situation where 
we have n dimensional state vector and p dimensional parame- 
ter vector, we define 


z 



Since k is constant. 



(17) 

(18) 


We adjoin Eq. (18) to Eq. (13) and use Eq. (17) to give 


i = g (t,z) ; z ± (o) = x ±Q , i= 1,2,.., n ; 

z i (o) = ? , i= n+1 , n+2 , . . , n+p 

The problem is to find the initial conditions z^o), i= n+1,.. n+p 
to minimize the least squares criterion 

R _ ni 

S = X- (z r - h (K>z (t ;z(o)))) Q r (£ r -h (t ,z (t r ;z(o)))) 
r=1 

where z (t;z(o)) is the solution of Eq. (13) corresponding to 
a particular set of z^(o), i= n+1,.., n+p. 

7. Certain short cut methods have also been presented in the 
literature (23,24)^ These methods are very useful for para- 
meter estimation as the estimates are obtained without solving 
the differential equations. In (23) ? the mathematical model 
is assumed to be linear in parameters and the model is comp- 
osed of M independent equations of the form: 



dC 1 

at" 
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+ *• + B k R 1 ,k< t > + ** B N R 1,N (t)?C 1 (o)=C o1 


U.U . 

at 2 = B 1 R j,1 (t) + •• + 3 k R j,k (t) +•• Vj,N (t)|C J (o)==C o2 (19) 


uo M 

dt = B 1 R M,1^ + ** +B k%,k^ ^ + *’ %%,N^ ;C M^ °) =C oM 

where the C's ( concentrations) are the dependent variables, 

the B*s are the kinetic coefficients, presumed constant, and 

the It's represent any desired function of C and t ( many of 

which may be zero). The Eq. (19) can be written in an abbre- 

dC T N 

viated form as B^Rj k (J= 1,2, . . , M) (19a) 

Eq. (19) has a unique solution if the right hand side of Eq. (19a) 

has a continuous uniformly bounded partial derivatives with 

respect to C. in the region of interest. Now the main difficu- 
3 

Ity in using Eq.(l9) directly to find the ( problem in the 
category 2 above) is that, as discussed above, the derivatives 
on the left hand side must be known and this cannot be done 
accurately. On the other hand the formal integration requires 
iterative calculations ( category 3). However, if each equation 
is directly integrated using the fact that the right hand side 
involves only integration of functions of the measured concent- 
rations as functions of time, good features of both categories 2 
and 3 are retained. 
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Thus, after integration from t Q to t^ Eq. (19a) 
may be written as 


o/n) - vv 


H 


2 — ) 
k=1 ^ 


E j,k « 


( i— 1 ,2 , . . ,p) 
U= 1,2,..,M) 


( 20 ) 


where P is the number of time intervals between the times, at 
which the concentrations were measured. Eq.(20) may be abb- 


N 

reviated as Y^ = C-y- C 1 j B k Xy^ (20a) 

k— 1 

The r.h.s. gives the predicted value of the modified dependent 
variable and observed values can be obtained from the experi- 
mental data. A least square criterion may be then used to 
estimate the unknown coefficients, B^. It is better to use 

sum of squares of the weighted deviations = (Yy-YylWy, 

/\ 

instead of sum of squares of the deviations (Y. --Y. .) , where ¥ . 

X J X J -I- J 

is the weight associated with each deviation. 


8. Sequential Experimental Design procedures for precise para- 
mater estimation in ordinary differential equations were first 
presented by Hosten and Emig (37)^ These designs are based on 
the ellipsoid region that surrounds the true values of the 
parameters. Box and Lucas^®^ suggested to select those exper- 
imental conditions which minimize the volume of the joint 
confidence region. This is realised when det V(b) or equi- 
valently, det (X^X) is maximized. This criterion is known 
as " Minimum Volume” criterion. V(b) is the variance-covariance 
matrix of the parameters. 
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Hosten 


( 39 ) 


advocated realising the greatest 


possible contraction of the longest principal axis of the 


hyper ell ip so id, rather than trying to attain the smallest 
possible dimen ion. It is done by performing those expert- 

•4 

ments which maximize the smallest eigenvalue of V (b) . 

This criterion which acts upon the shape of the joint confid- 
ence region is known as " shape " criterion. These designs 


determine the influence of the sampling range on the precision 
of the parameters and also that of the number of sampling 
points in each range. 



CHAP TEH 2 


P ARAMETER ESTIMATION IN ORDINARY DIEFEBENTIAL EQUATIONS 

The problem of parameter estimation in case of 
algebraic models has received considerable attention in the 
past and the results are quite well established on the basis 
of sound statistical principles. 

For a linear algebraic model the following basic 


formulas are well known 

Y = X| + § (21) 

b = (X T Z) _1 X% 

V(b) = (/X)" 1 ^- 2 (23) 

0:-b) T I T X (f -b) = s 2 pE ( p,n-p, 1 -«c ) (24) 


Eq. (21) relates the n observations Y of the dependent variable 

Y to the nxp values of the independent variables X Eq. (22) 

gives the least squares estimates b of the unknown parameters 

Eq. (23) and (24) are respectively the expressions for 

the covariance matrix of the estimates b and the expression 

for the joint confidence region E associated with b . Eq. (23) 

and (24) are based on the assumption that the experimental erro 

2 

rs are independent with zero mean and constant variance <c . In 
addition Eq. (24) requires that ^ is multinormally distri- 

p p 

buted i.e., £. N(0,I £ ) . s is the estimate of the error 

o 

variance 6 , obtained from the model. 
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An important problem in process analysis is the 
estimation of parameters in a differential model from experi- 
mental data. Often it is impossible to solve the differential 
equations explicitly and one is forced to estimate parameters 
directly from the differential equations. 

One way of going about this is to combine numeri 
cal integration of the differential equations with non linear 
least squares estimation of parameters, for which there are 
several algorithms available. 

Alternatively, the differential equations may 
be linearised and the parameters determined by a repeated 
application of linear least squares, an approach which forms 
the basis of the quasilinearisation technique as explained 
earlier. 

Variational principles may prove very helpful 
in circumventing the cumbersome job of repeatedly integrating 
a set of differential equations and function minimisation and, 
save a lot of computer time. Indeed, in many cases, when the 
equations are linear in parameters, these methods provide 
explicit estimates of parameters. Weighted residual methods 
fall in the category of these methods. 

Looking at the importance of parameter esti- 
mation in ordinary differential equations it was decided to 
use .different algorithms that are available for parameter 
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estimation in ordinary differential equations, and compare 
their suitability for a particular problem at hand. 

Pour different algorithms are compared for 
their fastness and reliability of the parameters estimated 
by them. The models, used for testing the efficiency of the 
algorithms represent batch reactor kinetic data for some of 
the most important industrial reactions viz, direct esteri- 
fication of terephthalic acid by ethylene oxide to produce 

f 2 6 2 7 1 

bis-hydroxy terephthalate and isomerisation of m- xylene ^ ’ 

In case of batch reactor, the models representing the kinetics 
of the reactions are given by first order or dinars 7 differential 
equations and the problem is essentially an initial value prob- 
lem. The choice of the algorithms was based on the different 
theoretical considerations that exist for solving differential 
equations and their ease of applicability. 

Runge-Kutta methods have an important place in 
the numerical integration of differential equations especially 
for initial value problems. An algorithm which uses a fourth 
order explicit Runge-Kutta-Gill algorithm was used to 
estimate the parameters of the models for two different kinetic 
schemes. These methods are quite sensitive to the stepsize of 
integration and the stiffness of the differential system. These 
problems are discussed, in detail, later. Because of its very 
frequent use nonlinear least squares analysis, was joined 
to the numerical integration scheme for the minimisation of 
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the sum of squares of the deviations of the predicted values 
of concentrations from the measured values. 

Another minimisation program based on the 
Box-complex method ^9) waS a q so tried to determine the 
effect of the efficiency minimisation on the estimated values 
of the rate constants. The algorithm was used in conjunction 
with a fourth order explicit integration scheme.. 

Orthogonal collocation methods w 7 are used quite 
frequently for the numerical solution of differential equations. 
The applicability of the technique. to initial value problems- 
e specially nonlinear first order differential equations leads 
to one point collocation method which is quite suitable for 
nonlinear coupled systems as the method is an implicit one. 

An algorithm based on an implicit second order scheme for nume- 
rical integration, coupled with nonlinear least squares analysis 
was used to observe an improvement in the values of the parame 
ters by using implicit equations. These methods are more stable 
and reliable at the expense of more computational time. 

The fourth algorithm used was based on varia- 
tional principles and used weighted residual method for para- 
meter estimation directly without solving the set of differen- 
tial equations. The method provides with an explicit estimates 
of parameters for models linear in parameters and even for non- 
linear models obviates the necessity of minimisation of an 
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objective function and numerical integration. A set of non- 
linear algebraic equations result which can be solved by using 
Newton-Raphson or successive substitution techniques. The method 
is extraordinarily fast and requires a negligible computer time. 

Let us define the problem first before suggesting 
the different ways to solve it. 

Consider the problem of estimating the vector of 
parameters © in the differential model 

dy 

** =£(£ ’ a) ( 25 ) 

with initial conditions Z = Z 0 > t = 0 

whe re 

(i) z is the n-dimensional state vector 

(ii) © is the m- dimensional vector of parameters 

(iii) f is the n-dimensional vector function of 
known model forms. 

(iv) t is the independent variable . 

As can be seen the equations are most general in nature and 
can represent continuous flow integral reactors with t repre- 
senting reactor length, Z. 

Suppose we are given a set of observations of 
Y(z) at p discrete values of t, assumed related to the state 

by 
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I r = I r + ~r ( r= 1,2,.. ,p) (26) 

where is an n-dimensional vector of observation errors. 

ALGORITHM 1 : 

Fourth order explicit Runge-Kutta-Gill method 
of integration use <j to solve the given set of differe- 

ntial equations and then the value of the state vector Y is 
used in Eq. (26) to calculate the vector of residuals. To obtain 
the estimates of parameters, the objective function given by- 
re si deual sum of squares is minimized by using a nonlinear 
least square minimisation routine. 

The R-K-G method used is described by 

dy 

ir e i 2 > 

assuming that 0 , the vector of parameters is known 
^n+1 = 2n + | ( K 1 +2K 2 +2K 5 +K 4 ) 

where K^ = f (Z n ) (27) 

K 2 = £ (Z n + \ ix K 1 ) 

K 3 = f (^ 1 +iK 1 +BK 2 ) 

R 4 = f (l n +CK 2 +DZ 3 ) 

where A,B,C and D are R-K-G constants given by 

A= J|=1 B= 2.-/I 0 =-=/T ^ Uf 

Once the vector of state variables at different time values 
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is known, residual sum of squares is evaluated which is a 
measure of discrepancy between the observed values and the 
model values is calcuated 



This objective function is minimised by using a nonlinear 
least squares minimization routine. 


In the flow diagram that illustrates the algori- 
thm: 

y rs - observations for the multiresponse system with 
r= 1,..,p and s = 1,.*,n 

f - vector of functions of v and © 

s _ 

9 0 - vector of initial guesses of parameters 

y - initial conditions for the n-responses 

SO 

t - discrete time values at which measurements are done 
r 

p - number of observations. 
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Fig. 2 Parameter estimation in ODEs by explicit BKG method. 
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Non Linear Least Squares Analysis : 


A modification of Gauss- Newton method as 
suggested by Grill and Murray is used in the present work for 


the solution of nonlinear least squares problem 


( 28 ) 


The 


algorithm is concerned with the location of a point x which 


minimizes the sum of squares of nonlinear functions 
F(x) = ^ (f. (x)) c , x £ S 


,n 


m 




(29) 


i=l '~l 

The gradient vector g (x) and Hessian matrix G (x) of F(x) 
are given by 2J(x) T f(x) and 2 ( J(x) T J(x) + B(x)) respectively 
where J(x) is the mxn jacobian matrix of f(x) whose ith row 

is Vf.(x) = ( af i^x 2 , B(x) 

m |< 

= X- f *( x ) G. ( x ) and G. ( x) is the Hessian matrix 
i=1 1 1 1 

of fj^x) (F(x) assumed to be twice-continuously diff erentiaHe) \ 
For uncons trained minimisation problems where the Hessian 
matrix of second derivatives can. be calculated, Newton' s Method 

\ , * ’ V 1 , 4 ■ l > 1 : ' | 

can be used. The method constructs a sequence of vectors 
£ x^ such that X (k+D = x (k) +oC (k) p (k) ( 30 ) 

where is a scalar steplength and p|^) f the direction ■ 

of search, satisfies the equation 

G(x (k) ) P^ k) = - g (x (k) ) (31) 

The special form of the Hessian matrix and gradient vector 
can be used in the Newton equation to give the equivalent 
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form (J(x (k) ) T j( x (k) ) + L(x (k) )) P^ k) = -J (x (k) ) T f(x (k) ) (32) 

The Gauss-Newton Method exploits the special structure of the 
Hessian matrix and gradient vector which occurs with least 
squares problem. The method computes the direction of search 
as the solution of 


J(x (k) ) T JU (k) ) p^ k) = -J (x^) 1 f (x (k) ) 


(33) 


These equations are obtained by neglecting the second deriva- 
tive matrix B(x^ k ^) in Eq.(32). The Gauss-Newton method is 
intended for problems where j B(x) jj is small compared to 
j| J(x)'*" J(x)j such as the so-called " small residual problem" 
where f(x) — t>~o as x — »x*. 


Eq. (33) are the so-called " Normal Equations" 
for the linear least-squares problem. 


minimiz e 
P 


J(s (k) ) p + f(*< k >) II 2 


( 34) 


that value of p is selected which minimizes the linear least 
squares problem and has least Euclidean length . Further insi- 
ght into the method can be obtained by ..ref erring to (28) . 

ALGORITHM 2 ; 

The method uses a fourth order Runge-Kutta- 
Gill scheme for numerical integration of differential equa- 
tions coupled with the complex method ( constrained simplex 

(29) 

method) proposed by Box v ‘ . The complex method searches 



30 


for the maximum value of a function f (x^j.. t x n ) subject to 
m constraints of the form x^. ^ h^, k= 1,.. , m, where 
x n+ i, * • t x m are functions of x^ , . . , x^ and the lower and 
upper constraints and h^. are either constants or functions 

( 30) 

of Jj,.‘ i x^ It has been developed from the simplex method w ; 

0 q 

It requires an initial point x^ , . x n which satisfies all 
the m constraints to start the search for the maximum (mini- 
mum with-f) of a nonlinear function. In this method k^s n+1 
points are used, of which one is the given initial point. The 
further k-1 points required to set up the .‘initial configu- 
ration are obtained one at a time by the use of pseudo-random 
numbers and the ranges for each of the independent variables, 
viz., Xj_= g^+r^Ch^-g^) where r^is a pseudo-random deviate 
rectangularly distributed over the interval (0,1). A point 
so selected must satisfy the explicit constraints, but need 
not satisfy all the implicit constraints. If an implicit 
constraint is violated the trial point is moved halfway tow- 
ards the centroid of those points already selected ( where 
the given initial point , is included) . Ultimately a satis- 
factory point will be found. Proceeding this way (k-1) points 
are found that satisfy all the constraints. The method uses 
cC, the reflection factor as 1.3 and k, the number of vertices 
for contrained problems, as 2n. 
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ALGORITHM 3: 


This algorithm is well suited for the solution 
of general first order initial value problems and can handle 
stiff systems as well. A brief description of the method will 
give a deeper insight in to the use of the method. 

let us first consider scalar differential equa- 


tion 


^2 = f ( y) y( x ) =y 

dx n' J n 

where f(y) does not depend explicitly on x. A Taylor series 
from (x n ,y n ) is 

y n+1 = y n + hf 4 h 2 fyf + £ h 3 (fyyf 2 +fy f+ ..) (35) 

where f = f(y n ) , f y = (|~) y ^ etc. 

How a one point Gauss method with collocation at x=x n +h/2 and 
extrapolation to y n+1 gives 

2(yi/ 2 -y n ) = & f (yi/ 2 ) (? 6 ) 

y n+ i = yC^+ii) = y n +2 (y-|/ 2 ” y n) (37) 

The method is 0(ir ) in y n+ i but one nonlinear equation 
2(yi/ 2 -y n ) = b. f(y l/2 ) should be solved for y 1 / 2 * If f ( y ) 

is linearised from y n , we obtain 

2 ( y l/2” y n) = h ^ f ^ y n^ + f y ^ y l/2 “ y n^ ^ 8 ) 

or -1 

y 1/2 _y n = (2 ~ h V h f 


( 39 ) 
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and finally 

y n+1 = y n + O" 2 f y h) ' 1 M (40) 

Expanding the above equation yields. 

y n+1 = y n + h f 4 h2 f y £ + I h3 f y f ( 41 ’ 

Continuing next to M coupled differential equations we obtain 
an analogous expression 

£-,=(!- 0.5 h f y ) _1 h f (£) (42) 

£ n+ 7 £n (43) 


It is a semi-implicit Runge-Kutta method and works very well 
even for stiff equations and represents one of the best methods 
to solve an initial value problem. To solve an extraordinary 
stiff differential equation an extension of the above method 
may be used which will have more favorable damping properties: 


K 1= ( 1-ahfy)" 1 hf 

(44) 

K 2 = (1-ahfy)” 1 K 1 

(45) 

y n+1 = y n + ^1^1 + ^2^2 

(46) 


where the constants a,R^ and k ave the values of 
a = 1- JT/2 R 1 = 1- JTs Eg = sfT 

ALGORITHM 4 : 



Weighted Residual Methods: 

A vector of equation residuals R is defined for 
the given set of differential equations by 
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R = 


d£ 

dt 


t <a>2) (47) 

the solution ofthe set of differential equations satisfies 
the condition R = 0 identically for all t ^ 0. It must also 
satisfy the condition 

S 

(48) 


J 


G R dt 


=0 


a 


where G is an nxn weighting matrix and t and t- are any t 0. 
The various weighted residual methods differ only in their 
choice of G* 


The Residual Least , . Squares Method (RL£H) : 


mm 

6 


Instead of Eq. (48) we may use 
! = J b R T | E dt jr 


(49) 


where ^ is an nxn weighting matrix. Minimisation of S leads 
to m normal equations: 


<5 S 

ae. 


G= 


= 0 

(i= 1,2,.. 

> m ) 

choose — W = I , 

then 

where 



^ff 

* f 2 

^ f n 

a© 1 


d 

df 

df 2 

f 

n 

d© 2 

d© 2 

*• 

’ 

a©J *•** 

d!i_ 

d e m 

b± 2 

TFQ~ **•* 

m 

* 9 m 


(50) 
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combining Eqs« (47)& (48) we obtain 



& d Z 


J £ £ (l»£ ) dt 

t 


a 

substituting observations *Y for the state £ Eq* (51) 
to a system of nonlinear algebraic equations in the ©^ 


(51) 


reduces 


T lr A Zr = I (l r > 2 ) &-*r 


(52) 


However, if the function f is given by f = A (^) . © 

then Eq. (52) becomes explicit in © : 

P*~ 1 m _a P~ 1 m 

£ = ( r (A 1 A) A.t r ) 1 . ( A T A £r ) (53) 

r=o r=o 

Approximate Analytical Solutions: 

The differential equations of first order can 
be solved analytically but for a set of ordinary differential 
equations that are coupled, only approximate analytical soluti 
ons are possible. These solutions can be utilised very well 
to prepare a general purpose program for parameter estimation 
in the models governed by linear differential equations. 

A model governed by first order linear differen- 
tial equations is given by: 

; J (o) = I 0 ( Given) 
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This vector equation can be solved analytically 
A( t-t 0 ) 

to give Y (t) = e Y^ = exp (A (t-t Q )) Y Q (54) 


The exponential matrix is defined by the infinite series 

^ I + At + A 2 t 2 A 5 t 5 

e = ~ 2 \ + ~ 3 ! + * * ' 

let us define the state transition matrix as 

A (t-tj 

| ( t-t Q ) = e 


(55) 


(56) 


Therefore, we have 


I (t) = g (t-t Q ) Y o (57) 

where we may use a truncated series expansion of the expo- 
nential matrix as an approximation to ^ matrix. This matrix 
is made up of terms containing parameters and thus is a 
" constant" matrix. For different values of t, one gets the 
state of the system which undergoes transition . 


An algorithm would make the implementation 
of the procedure more clear (Fig.3)* 

The responses were . calculated both based on 

(t-t Q ) where t is constantly increasing and also based on 

( t . - t . ) where %.+1 and% are consecutive two time values. 

1 4-1 1 

The second procedure gave better estimates as was evident 
from the sum of squares at the end of the program. 

The values of the parameters obtained by using 
this analytical procedure may be used with better degree of 
confidence than those obtained by using purely numerical 
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Pig. 3. Algorithm for the implementation of analytical 


solutions 
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techniques. A still "better approximation to the optimum 
may be obtained by using a different objective function 
which would take into account the covariance of the residuals, 
especially in the case of multiresponse models. 

Formulation of the Objective Function : 


Another objective function based on the discri- 

( ?LA \ 

mination criterion of Box w ' was constructed and used in 
parameter estimation. The results were astoundingly better 
as was evident from the comparison of residual sum of squares 
with that of this objective function. 


From the theory of linear regression analysis , 

T 

it is known that the determinant of the matrix X X, where X 
is the design matrix of independent variables as pointed out 
earlier, can be used as a measure of the reciprocal of the 
variance-covariance matrix of the parameters in case of models 
governed by Y = X B . Thus, the determinant of ( X T X) ^ may 
be used as a measure of the deviation of parameters from their 
true values and thus, may serve as a very relevant objective 
function. The objective function used is given by: 


£ (yn-yil) 2 Z_ 1 ( y il” y il)( y i2” y i2) 

H / A v ✓ A v , A v 

£ , < yi2- y i2> ( y n - y n > f =1 ( yiz-yiz) 


• • ^ y A1~ y i1 ^ y iH y iM^ 

2 U A A 

* * * F-l ^ y i2_y 12 ^ y iM“ y iM^ 


N ^ ^ 

? ^ y ii~ y ii ^ 

1=1 


N A 

•^_1^yiM -y iM^ 
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The algorithms discussed above along with the 
objective function constructed above were used on the real 
time data available from the literature ^ 26,27) so the 

results obtained may be of some practical significance. The 
data used is reported. 


EXAMP LE 1 : 

Direct Esterification of Terephthalic acid with 
Ethylene Oxide : 

This reaction is homogeneously catalysed by 
triethyl amine (atertiary amine) and the reaction mechanism 
is given below. 


The reaction between ethylene oxide and tertiary 
amine gives rise to the formation of intermediate (quarternary 
base salt) which is immediately neutralised by Terephthalic 
acid. 

0 ^ + 

Et^N + CH 2 CH 2 fast Etj ITCH 2 CH 2 0(-) 

tri ethyl ethylene 

amine oxide 


+ 

Et 5 NGH 2 CH 2 0’"+ 


CooH 



COOH 


+ C oo~ 

fast Et 3 lTCS 2 CH 2 OH + 

• COOH 



Terephthalate anion thus formed reacts with ethylene oxide to 
give 2 hydroxy e\hyl terephthalate. 
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Coo 0 

(O) + CH 2 CH, 

coon 

^OOC«,C HjO~ 

o 


X- 


H 


+ 


WOH 



+ CH 2 

^ooch^ch^o" 

Hydroxy ethyl 

terephthalate 

ion 


-CH, 


slow 


fast 


slow 


CDOCW.COO 


© 


CO OH 


COcCH^X'ti.Ot! 

o 


Coo« 

^^OOC H jCHjOH 

0 

COOC^CHjO 


CoOCvt 5 t H-> OH 


H 


o 

Coo 


>fooc H^c H^OH 

fol 

COOCH^C Hj_OH 

8He.*r 


Purther, the equilibrium constant of the salt 

formation from a weak acid and a weak base is expressed as ; 

k, 

RCOOH + R‘1TH 2 — j — RCO~ + R’NH+ 


k . J R000 ) (R'NE 3 ) ) 
(RCOOH) (R*HH 2 ) 


k 1 


The above equilibrium relationship holds good in the absence 
of solvent Butyl alcohol that is used. In the presence of 
n-Butanol the equilibrium relation is expressed as k.= 


where kg is the ionic product of the solvent, k^ & kg are 
the dissociation constants of acid and base respectively. 
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The rate expressions governing the mechanism 

discussed are: 

-d(A)/dt = 2k 3 (A)(B)(F)/(2(A)+(C)) (59) 

-d(B)/dt = (B)(1) (2k 3 (A) +k 4 ( C) )/( 2( A) +(C) ) (60) 

d(C)/dt = (B)(1) (2k 3 (A)- k 4 (G))/(2(A)+(C)) (6t) 

d(D)/dt = (B)(1) (k 4 (C) /2 (a) +(C>) (62) 

where (A) = concentration of terephthalic acid in gmol/l 

(B) = concentration of ethylene oxide in gmol/l 

(C) = concentration of hydroxyethyl terephthalate 

in gmol/l 

(D) = concentration of bis-hydroxy terephthalate 

in gmol/l 

These four equations represent reactant and product distribution 
as a function of time. 

These equations are solved for the optimum values 
of the parameters k 3 and k 4 which are the reaction velocity 
constants for the formation of HET and BHET respectively. 

EXAMPLE 2 

Isomerisation of m- Xylen e: 

The reaction scheme for the isomerisation can be 
represented as 

(i) Series Reaction : p-Xylene * m-Xylene 

- — o-Xylene 

(ii) Coupled Reaction: o-Xylene - — »■ p-Xylene - — ^ m-Xylene 
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1 




3 


p m o 

k| , k^,k^,k^ are reaction rate constants. 

The differential model representing this reaction mechanism 
would bet 


dC 1 
d i 

-k'sC 1 + k^ 

SC 2 


(63) 

dC 2 

dt“ = 

k| SC 1 - k' 2 

sc 2 - 

^sc 2 + k| sc 3 

(64) 

dC 3 
di“ = 

k^ SC 2 - k| 

SC 3 


(65) 


where at any time t, C-j , C 2 , are the concentrations of 
p, m, and o-Xylenes respectively and S ia the active surface 
area of the catalyst. The rate constants icj^s are based on 
unit effective surface area. For the same catalyst sample 
used in different reactions the value of S remains unchanged 
and so it can be coupled with the rate constants k£s to 
yield the rate equations .as follows: 


dC 


1 


dt 

dC 2 

“dt 


= -k 1 c 1 + k 2 c 2 


- k 1 C 1 k 2 C 2 -k^C 2 


dC, 

dT 


where k = k'S hr . 


k^C 2 — 

-1 


( 66 ) 

(67) 

( 68 ) 
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The equations are also solved by the methods 
explained and discussed in the earlier chapter for the 
optimum value of the four reaction velocity constants 


kj_: i= 1 ,2,3 ,4. 


The kinetic data used for the estimation of 

/ 2 ^ 2 r 7 \ 

parameters in two different models elucidated above are: * 




43 


MODEL 1 ( from Beferen.ce 26) 

TABLE 2 

experimental conditions 


Run 

No. 

wt. of 

catalyst 

(grams) 

Cone. of 

catalyst 

(gmol/1) 

Temperature 
of the run 
°C 

Initial 
TP A 

(grams) 

Total 

volume 

(ml) 

1 

9.9974 

0.185 

60 

44 

535 

2 

2.7030 

0.043 

60 

1 66 

624 

6 

2.7159 

0.0434 

60 

166 

620 

7 

3.7990 

0.061 

ioo 

166 

620 

8 

9.2930 

0.172 

60 

44 

535 

10 

10.1699 

0.160 

60 

166 

630 
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TABLE 5 

EffiEMM DATA 



Sample Time 
No. (min) 

TP A 

(gmol/l) 

Bo 

( gmol/l) 

nsr ~"” 

(gmol/l) 

BHET 

(gmol/l) 

Run 

1. 

0.0 

0.476 

1.750 

0.000 

0.000 

No. 1 


2. 

15.0 

0.416 

1 .670 

0.027 

0.03 5 


3. 

30.0 

0.393 

1.580 

0.034 

0.049 


4. 

60.0 

0.333 

1.373 

0.051 

0.092 


5. 

90.0 

0.239 

1.363 

0.071 

0.165 


6. 

120.0 

0.186 

1 .257 

0.079 

0,211 


7. 

180.0 

0.156 

1.000 

0.099 

0.220 

Run 

1. 

0.0 

1.590 

1.545 

0.000 

0.000 

No. 2 


2. 

15.0 

1.570 

1.543 

0.007 

0.006 


3. 

30.0 

1 . 540 

1.518 

0.010 

0.037 


4. 

60.0 

1.530 

1.395 

0.010 

0.049 


5. 

90.0 

1 .510 

1.318 

0.010 

0.056 


6. 

120.0 

1.470 

1.300 

0.012 

0.068 


7. 

180.0 

1.460 

1 .268 

0.017 

0.074 

Run 

1 . 

0.0 

1 . 600 

3.140 

0.000 

0.000 

No. C 


2. 

25.0 

1.430 

2.800 

0.010 

0.163 


3. 

40.0 

1.365 

2.680 

0.020 

0.234 


4. 

70.0 

1.333 

2.630 

0.017 

0.257 


5. 

100.0 

1.240 

2.440 

0.015 

0.348' 


6. 

140.0 

1.120 

2.240 

0.020 

0.461 


7. 

180.0 

1.017 

2.050 

0.022 

0.561 


Contd 
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Sample Time 
No. (min) 

TP A 

(gmol/l) 

E0 

(gmol/l) 

HET 

( gmol/l) 

BHET 

(gmol/l) 

Run i . o.O 

1.560 

1 .910 

0.000 

0.000 

No. 7 

2. 30.0 

1 .510 

1 .835 

0.01 5 

0.031 

3. 55.0 

1.340 

1 . 510 

0.037 

0.183 

4. 70.0 

1.111 

1 .052 

0.049 

0.400 

5. 100.0 

0.921 

0.754 

0.054 

0.585 

6. 130.0 

0.689 

0.302 

0.064 

0.807 

7. 180.0 

0.572 

0.010 

0.071 

0.917 

Run 1 . 0.0 

0.462 

1.560 

0.000 

0.000 

No. 8 

2. 15.0 

0.419 

1.493 

0.013 

0.030 

o 

• 

o 

K\ 

• 

K\ 

0.400 

1.409 

0.015 

0.048 

4. 60.0 

0.360 

1.395 

0.017 

0.085 

5. 90.0 

0.333 

1.364 

0.022 

0.107 

6. 120.0 

0.246 

1.323 

0.049 

0.167 

7. 180.0 

0.196 

1.168 

0.051 

0.215 

Run 1 . 0.0 

1 .500 

1.560 

0.000 

0.000 

No. 10 . 

2. 25.0 

1 .470 

1 . 510 

0.005 

0.024 

3. 40.0 

1 .450 

1 .448 

0.011 

0.042 

4. 65.0 

1 .410 

1 .395 

0.017 

0.072 

5. 100.0 

1.272 

1.240 

0.019 

0.209 

6. 140.0 

1.166 

0.900 

0.011 

0.323 

7. 180.0 

1.106 

0.800 

0.017 

0.377 
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MODEL 2 ( from Reference 27) 

■TABLE 4 

. EXPERIMENTAL CONDITI ONS 


Data 

Set 

No. 

Catalyst Reaction 
composi- Temperat- 
tion ure (°K) 

Catalyst 

loading 

(gm/ml) 

Catalyst Stirrer 
particle speed 
size (mm) (rpm) 


1 1 

.5$ Ni-HM 578 

0.02 

1.6 ' 

250 



2 1 

.5$ Ni-HM 533 

0.02 

1.6 

250 





TABLE 

5 







EXPERIMENTAL DATA 





Sample Time 

p- Xylene 

m- Xylene 

o- Xylene 

Toluene 

TMB 


No. (hrs. 

) gmol/l 

gmol/l 

gmol/l 

gmol/l 

gmol/l 

Run 

1 

0 

0.0000 

4.0000 

0.0000 

0.0000 

0.0000 

No. 1 









2 

1 

0.1252 

3.7512 

0,1068 

0.0072 

0.0080 


3 

2 

0.2172 

3.5700 

0.1832 

0.0144 

0.01 52 


4 

4 

0*3960 

3.2080 

0.3352 

0.0304 

0.0304 


5 

6 

0.5433 

2.9068 

0.4624 

0.0448 

0.0420 


6 

8 

0.6700 

2.6448 ' 

0.572 

0.0560 

0.0572 


7 

9 

0.7200 

2.4988 

0.620 

0.0808 

0.0804 


8 

11 

0.7630 

2.3580 

0.704 

0.0840 

a. 0860 


9 

14 

0.7736 

2.3320 

0.7144 

0.0900 

0.0900 


10 

16 

0.8520 

2.2272 

0.7680 

0.0856 

0.0872 


Contd. . . . 





1 


0.0000 


Hun 1 
No. 2 


2 

2 

0*0380 

3 

4 

0.0704 

4 

6 

0.11 60 

5 * 

8 

0.1650 

6 

9 

0.2360 

?■ 

11 

0.2600 

8 

14 

0.3280 

9 

16 

0. 3504 


4*000 

0.0000 

0.000 

3*938 

0.0296 

0.000 

3.870 

0.0600 

0.000 

3.762 

0.1000 

0.0108 

3. 668 

0.1420 

0.0128 

3.532 

0.1930 

0.0196 

3.552 

C .3080 

0.0220 

3.314 

0.2980 

0.0248 

3.287 

0.3060 

0.0280 


iOOO 

0.000 

0.000 

0.0112 

0.012 

0.0196 

0.0216 

0.0268 

0.0280 
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CHAPTER 3 

CONVERGENCE AH ID STABILITY OE PI EFBRENT I AL ECU AT I OH S 
5 a . INITIAL ESTIMATES OE PARAMETERS : 

Most of the optimisation routines require a 
good initial guess of parameters to attain the optimum value 
of the objective function. The solution may overshoot the 
optimum or convergence becomes excessively slower if an 
initial guess is quite far from the optimum. This problem 
would be much more severe in case of differential models 
where divergence due to poor initial estimates may get 
compounded with the instability and divergence due to the 
nonlinearity of the equations and poor step size used in 
integration. 

A variation of linear least squares 
technique is used for differential models for multiresponse 
systems to provide near optimal initial estimates of the 
parameters involved in the model: 

In case of linear algebraic models 

Y = I £ + 8 

The best linear unbiased estimate ( BLUE) of parameters is 
given by 

b = (A)" 1 A 


( 69 ) 
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where Y is the vector of responses at different times Z is 
the constant matrix of independent variables and b is the 
BLUE of § . 


In case of differential models we have: 


dY 


dt 

Y 


= Z 


= Y + € 


(70) 


The vector Y. at different time values is observed. 


dY 

In this model, — 

dt 


is replaced by n.Y 

It" 


-which 


is nothing but the slope of Y v/s t curve at the average 

value of two consecutive concentrations at different times. 

The design matrix Z - as it is known as - is calculated for 

each time interval at the average concentrations (keeping' 

in view the mean value theorem of calculus ) . The BHS of 

T I 

Eq. (69) is evaluated in two segments : Z Z and Z Y and 

summed up over all the time intervals to form two different 

sum matrices respectively , The former sum matrix is 

inverted and multiplied with the latter to give a very 

closely optimal initial values of parameters: 

P“1 m 1 P-1 m 

b=(H II)" (X t:± L) (7 1 ) 

r=o r=o At 

Where p observations are made of different reactant and 
product concentrations involved. A program is attached in 
the appendix for carrying out these calculations. 
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3t>* STABILITY GHABACTERISTICS OF DIFFERENTIAL MODELS ; 

The integration methods used to solve a set 
of differential equations are susceptible to instability which 
is largely caused by a large stepsize used and, sometimes, 
due to highly nonlinear nature of the equations. 

We will try to understand the phenomenon 
caused due to large step size. Let us take the test equation 

|§ = -Ay y(o) = 1 (72) 

where )s. is real and positive. Let us, then, write the 
solution as the sum of the exact solution y and an error € . 
We put this equation into Bq. (72) and note that the exact 
solution satisfies the differential equation , too. Then 
the error satisfies; 

Sf = (73) 

We examine the error in successive time steps by looking at 
■^ n = £(t n ) and £- n+ i • An integration method is stable if the 
error decays in successive time steps. Because of round off 
errors the computer never solves equations exactly. If the 
scheme is unstable this round off error grows with successive 
time steps and soon swamps the solution. 

We have used fourth order Ilunge-Ku 1 1 a- G-ill 
method and its variations to solve the system of nonlinear 
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first order initial value differential equations* So, we will 
calculate the maximum step size that can be used to ensure 
stability in the integration of our system. 

The R-K-G method is expressed by the algorithm 

K i = h f (y ,t ) 

1 w n* n' 

K 2 = h f (t n +l/2 h, y n +l/2 K 1 ) 

( 74 ) 

= h f (t n +l/2 h, y n +AE 1+ BK 2 ) 

K 4 = h f (t n +h, y n +CE 2 +IK 3 ) 

y n +i = y n + i/6 (e 1 +e 4 )+ 2/6 (bk 2 +h: 5 ) 

As the analysis in case of non linear functions is carried 
out by linearising them around y n ,t n , we will carry out the 
analysis assuming that equations are represented linearly in y. 

ti = ->- y 

Then, K 1 = -Xhy^ 

K 2 = - \h (y n -1/2 X hy n ) 

= -Xh (y n +(“ if-) ( - Xhy n ) +( ■ 

(- Xh (y n -l/2 Xhy n )) 

K 4 = - Xh (y n - J f (- Ah{y n 4 \Hy n ))+(lW| ) Kj ) 

After simplification, we get : 

K, = - Xhy n 
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K 2 = - ^hy n ( 1-1/2 Xh) 

Kj = - Xhy Q (2- Xh + (^ji- ) X 2 h 2 ) 

k 4 - - Xhy n (i +( J!=!^I)Xh + A 2 4xV) 

Therefore, substituting in 

y n+1 = y n + h ( K 1 +K 4) + § (BK 2 +DK 3 ) (77) 

gives 

y n+ , = ((.-a- /2) m + ( 

- £ 4 h 4 ) (78) 


Let us name the second terms on R.H.S of Eq. (78) as stability 

term 

then# 


as the 


y n+1 = y n ( 1+ 1 ( stability term)) 
exact solution satisfies the differential 
^n+1 = ^i ^ 1+ i (stability term)) 


(79) 

equations, 

(80) 


Bow, for a method to be stable, we should have 
-=> | ( Instability term) )1^ 1 


6 , 
n+1 

£ nT 


- 12<C" stability term 4^0 

Using these limits on the stability term, the term 
evaluated which turned out to be 


4.1 (si) 
(82) 
( 83 ) 

h was 


\Xh\ ^ 2,8 


(84) 
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where X is the maximum eigen value of the jacobian matrix 
in the equation 


ay- 


i 

dt 


= f i ^ (y -y ) 

3 ay, y 3 y jn ; 

J 


(85) 


ay 


_ A , the jacobian . 

3 c . A ROTE OK STIMES3 OF DlfTERENTIAl EQUATIONS^ 2 ^ : 


As we have seen in case of differential equations 

governed by ^ % 

"W = 4 Z z(o) = G (given ) 

4 is the matrix of coefficients 

That, the largest step size in governed by A 1 ( where p has 

' A maf 

different values depending upon the method and X ma3; is the 
largest of the absolute magnitudes of eigenvalues. 


We have, the unfortunate situation with system 
of equations that the largest stepsize is governed by the 
largest eigenvalue and the final time is unusually governed 
by the smallest eigen value. 


We define the stiffness ratio, SR as 

max \ Re X ^ 

_ i 

~ min l Re \ t 
i ’ i 


( 86 ) 


Typically SR = 20 is not stiff, SR = 10^ is stiff and SR= 10 
is extremely stiff. 




CHAPTER 4 


BELT ABILITY ANALYSIS 

In comparing various algorithms for parameter 
estimation reliability analysis can give us the degree of 
confidence that can be attached to the parameter estimates. 

A procedure of computing confidence regions or 
confidence interval s-e specially suitable for differential 

( 33 ) 

models- is presented as described by Rosenbrock and Storey v . 
Analytical Developmen t: 

Consider two adjacent solutions X and X + of 
= f ± (X ,a,t) i= 1,2,..;n (87) 

where a represents a p-vector of parameters. 

Starting at t=0 from X (o) = C , J|_( o) =0 where the perturba- 
tions are due to small changes =6 in a*: 

|sf (| ,a*,t) ; X(o) = C (88) 

X = f ( X+£ « a*+ oi ,t) ; X(o) =C 4 (o)=0 (89) 
Expanding Eq. (89) in a Taylor series, and keeping only terms 
of the first order inland oi . , we obtain 

X + Jr. = f(X,a*,t) + A I. + B < 

-i = a 4 +b * ; k (°) =o 


(89a) 
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where £ = (A^-) = §—■ f ± (2, a* ,t) 

I = ( B ij } = ir. f i u »a*» t ) 

3 

the matrix A is nxn and B is nxp; both matrices are func- 
tions of t but not of or . 

Now, consider the matrix differential equation 

Y = - A’ I; I(o) = I 

where Y is an nxn matrix and I is the identity matrix. We have 

It r & = r 4 + 1 ' i 


=- Y'4| +1' ( | |_ + B £) = Y'B 

cancellation of terms containing 1| is the motive for intro- 
ducing the adjoint equations. Then since_§ ( o) = 0 

Y'(t r ) & (t r ) = j r |*<.t ) |(t) <dt 

° (r=1,2,..,n) 

= 5 (t r ) £ 

where nxp matrix 0(t r ) is C(t r ) = j 11, Y'(t) B(t) dt 


The matrix Y*(t ) is assumed to be nonsingular for all t r . 


and then we have:|s (t^J = (Y'(t r )) C(t r ) 5 ("^ r ) — Sa Y 

where D is nxp. If follows that 


=r 


C D^(t )) = i (a*,t )) 


'ij'/r'V S»aj x 
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Thus* gives the dependence of X(t ) upon variations 
of a around a* which is what we set out to find. 

Consider next the function F(a*+ ) as given by 

R ** 

F(a) = ^ «‘(X(a,t r )^ W r |Y(t r ) -g (X(a,t r ) )j» 

Let c < he chosen to make F( a*+ a minimum, so that a*+ =a. 

Then for i= 1,2,.,, p. 

R r 

II {(r(t r )- g* (X(a* + £ , t r )) ¥ r (Y(t r )-g 

(X(a*+ 4 . »t r ))^ =0 


or 


R 


-JZ7 («r + V r“ ^D'spt^gr+qr-Sr-Vr^^ 

where g y = g (X(a*,t r )) 

gr •'% V» 


=0 


■ D ' x ' o 

and D is given by Eq. ( 


3V 


The matrix G^ is mxn, D r is nxp. 


Thus, we have 


R 

y =r = =r ' ' 7 r =r=r 


2’ S W„ ( Hi - SB «t ) =0 


r=1 

(£ SJ I r g r Br) Sf = II Br 2r Sr > 

r=1 r=1 

or R 

i:4=x g. a ir 'v 

r=1 - 


( 90 ) 
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If the symmetric matrix H is nonsingular above 
EQl* determines from the . The expectation of vC will 
be zero, and we are interested in the covariance of the 
elements of <K \ that is, in the expected value Pof .<./ 1 . 
By Eq. (90) 


R R 

E(g ^'1 ) = E(r n 

r=1r=1 


D' G' W n n W* G D ) 

—TV — -T* — T» -r — -T* ‘ 


Whence by using E (n n* ) = 0 r ^s and E (n n* ) 

X S XT 


M 

=r 


we have 


EPH 


R 

r^1 


D* G* W M W G ' D 
=r =r=r=r =r Sr Sr 


This equation, determines P if H is nonsingular . By Eq. (90) 
oC is a linear function of n r so that its probability density 
will be Gaussian and will be given by: 


1 

( 2 ' 7') p / 2 ~/ rg \ 




If now b is a unit p~vector the linear function b'y., of^C 
will have Gaussan distribution with variance 

= V P b 

The confidence interval r (b) associated with b 1 a is then 

r(b) =K<T- d (92) 


All the equations used sofar can be put into a 
computationally convenient form 
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X 

a 


~ -^(2» I ^(o) = C 

R 


(93) 


min j il* (t ) - g « ( X( a, Z (t ) 

a f =1 0 r J r t r (94) 

- £ (X(a, t r ))| 


4 = 

^ij'* . ^i a ) 

3 

(95) 

B = 

(B, ,) = (-fL. % (X,a* t)) 

•*- tj c* 

(96) 

¥ 

D = 

A(*) g + 5 (t) ; D(o) =0 

(97) 

8 W 

i! 

s (v 

(98) 

ir= 

g(t r ) = ( 8 13 (t r )) = f- gl (I(i, t r ))) 

R 3 

il, Sr 2r ir ' 1 2r Sr 

(99) 

H = 

( 100 ) 

P = 

B ' 1 

( 101 ) 

2 

b* P b 


b ~ 

( 102 ) 


& v, = 


In applying these formulae » it is assumed that 
a is found by using standard optimisation methods. This 
will involve many successive integrations of Bq. ( 93 ) with 
different values of a . 


When a has been found it becomes possible to 

estimate df, and so to find the confidence Intervals r(b), 
b ~ 

For this purpose a further integration of Eq. (93) is begun, 
using a = a . At each step of integration 4 and B are 
evaluated: Thus, the solutions of Bqs. (93) and (97) is 

carried out simultaneously, and D is therefore available 
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whenever it is required. Since X is available G can also 
be computed whenever it is needed. 

At each time t r , when an observation is available, 
the value of JQ^, fi-j, M r ^ £ r can be computed and added 
to the running sum which forms H. Once H matrix is obtained, 
by inversion P can be obtained or analysis can be carried 
out by finding eigenvalues and eigenvectors of H ; Eigenvalues 
of P are reciprocal of eigenvalues of H and eigenvectors are 
same. These eigenvalues and eigenvectors reveal, the situation 
completely. 

4a. EXAMINATION OP RESIDUALS : 

Residual analysis is useful and valid not only for 
the linear regression models but also for nonlinear regression 
models and other forms of models. 

The residuals are defined as the n differences 

A, A 

e. = Y.-Y., i= 1,2..,n . Where Y. is an observation and Y. 
ill' i i 

is the corresponding fitted value obtained by using proposed 
model equations. 

Now while performing the regression analysis we 
have made certain assumptions about the errors; the usual 
assumptions are that the errors are independent, have zero 
mean, a constant variance £ 2 , and follow a normal distribution. 
Thus, if our fitted model is correct, the residuals should exibit 
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tendencies that tend to confirm the assumptions that we have 
made, or atleast, should not exibit a denial of the assumptions. 
After the examination of residuals we can conclude that either 
1 * The assumptions appear to be violated ( in a way that 
can be specified ) > or 

2. The assumptions do not appear to be violated. The ways 
to examine residuals are graphical as well as numerical and 
are easy to do; they are very revealing when the assumptions 
are violated. The principal ways of plotting the residuals 
e^ are 

1 . Overall 

2. In time sequence; if the order is known. 

3. Against the fitted values 

Statistics for examination of residuals : 


The plots we have discussed above are visual tests 
to check some of the basic assumptions underlying the estima- 
tion of parameters. We may, however, use statistics to 


do the same analysis: 

n 


T. 


PI 


T 21 ~ 


q 

X- i 

i=1 

n 2 a 

X °i X 

i=1 


(103) 

provides a measure of (104) 


the defect of changing variance in time sequence 
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11 ~ 

n 

X 

e i y i 

This should always be zero (105) 

i=1 




n 

/v 2 


12 = 

y 

p. Y. 

1 

provides a measure of defect 


®L=1 




due to the absence of linear and quadratic terms of time 
in the model. 



CHAPTER 5 


RESULTS AMD DISCIJSSIORS : 

A comparison of four different algorithms to 
estimate parameters in models governed by ordinary differential 
equations is carried out so as to help the future investigators 
in selecting a method for parameter estimation in kinetic 
studies involving dynamic responses. 

NUMERICAL PROCEDURE: 


The four algorithms discussed in chapter 2 
were applied on the experimental data presented in Tables 2-4. ■ 

In the first algorithm we require initial condi- 
tions to start the integration of the equations. Initial values 
of the parameters are also required. The step by step procedure 
is given below: 

( 1 ) Model Used: 

(a) d(A)/dt = - 2 k 3 (A)(B)(F)/ (2 (a)+(C)) 

(b) d(B)/dt = - (B)(F) (2k 3 (A)+K 4 (C)/(2(A)+(C)> 

(c) d(C)/dt = - (B)(F)(2k 3 (A)-K 4 (C))/(2(A)+(C)) 

(d) d( D) / dt = - (B)(F)(K 4 (C)/2(A) +(C)) 

where the symbols have the meaning as given in chapter 2, 
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Run No.1 


Initial conditions are Initial values of the 

lgmol/1; parameters (mol/l)“1 J_ 

min 


(A) 

= 

0.476 

(B) 

= 

1.750 

(0) 

= 

0.000 

(D) 

= 

0.000 

(P) 

= 

0.185 

Run Ro.2 

(A) 

= 

1 .590 

(B) 

= 

1 .545 

(c) 

= 

0.000 

(D) 

= 

0.000 

(P) 

= 

O.C43 

Run 11 o. 6 

(A) 

= 

1.600 

(B) 

= 

3.140 

(c) 

r = 

0.000 

(D) 

== 

0.000 

(P) 

= 

0.0434 


K 3 = 0.0107 
K 4 = 0.0497 


Constant for the entire run) 


K 3 = 0.0092 

K 4 = 0.0762 


= 0.0066 

K 4 = 0.1170 



Initial conditions are 
( gmol/l) 


(A) =1,560 

(B) = 1.910 
'(C> = 0.000 

(D) = 0.000 
(P) = 0.061 

Run Ho. 8 

(A) = 0.462 

(B) = 1.560 

(C) = 0.000 

(D) = 0.000 
(P) =0.172 

Run N o . 1 0 

(A) = 1.500 

(B) = 1.560 

(C) = 0.000 

(D) = 0.000 
(P) =0.160 


Initial values of the 
parameters (mol/l)~1 1_ 

min 

= 0.1310 
K 4 = 3.2240 


K 3 = 0,0072 

K 4 = 0.1670 


K 3 = 0.0102 
K 4 = 0.0720 


The R.H.8. of Eqns. (a) - ( d) is evaluated at 
the given conditions with the parameter values as mentioned 
above. This function is used as f (T) in the algorithm and 
the R-K-0 parameters E^K^E^ and E 4 .are evaluated. Using 



these constants, we can get the concentrations at the n«ext 
time value* The time interval between two consecutive obser- 
vations is used as the stepsize in the solution of the equations 
Once the concentrations of A»B,C and D> as predicted by the 
model are available, these are subtracted from, the observed 
values reported in the data tables* These residuals axe 
then squared and added to a running sum to give the value of 
the objective function for the minimization program. Chen the 
minimization program carries out search for the minimum and 
the parameter values are changed. Again starting from the 
same initial conditions, the above procedure is repeated to 
give . - the next measure- of the objective function which is 
again checked for the minimum. This way the search continues 
for the minimum objective function and once it is reached, the 
corresponding values of the parameters provide the optinaal 
estimates* A prograh listing is included in the appendix. 

In algorithm 2 , a very similar procedure is 
employed till the calculation of the residuals, then tine vari- 
ance— covariance matrix of these residuals is constructed. The 
procedure is discussed in chapter 2. Determinant of this 
matrix is found out by employing complete pivoting technique 
as given in the program at the end . This determinant serves 
as the objective function and is fed into the complex niaiimi- 
zation routine. The parameter values are updated and again 
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starting from the initial conditions, the objective function 
is evaluated till a local optimum is reached by the program. 

The values of the parameters corresponding to this situation 
are reported in the tables 6. 

In algorithm 3 , the equations are solved implicitly. 
A jacobian matrix J ac is constructed first. This is the 
matrix of derivatives of the functions on the R.H.S. of the 
equations. This matrix is multiplied by o.5h which is half 
the difference in the two successive time values. Now this 
matrix is evaluated at the initial conditions given above and 

then subtracted from an identity matrix of the same order, I 

\ 

as discussed in the chapter 2. The resulting matrix is inver- 
ted and multiplied by the matrix F (Y,0) which is obtained by 
substituting initial conditions and initial parameter values. 
This, is then multiplied by h again. The result is vector 
which is added to the initial conditions to give the concentra- 
tions A,B,C and D at the next time value. This way the calcul- 
ations are performed for all the time intervals starting from 
0 to 180 minutes. These predicted values are subtracted from 
the observed values and then the residual sum of squares 
evaluated. This is fed into the nonlinear least square s 
program as objective function. The resulting parameter values 
are then used along with the same initial conditions and the 
procedure repeated until optimum is reached. The optimum 
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is reachedi 

parameter values are reported. A detailed program is 
attached in the appendix. 

In algorithm 4, matrix A is first evaluated. 
Tliis matrix is has concentration terms as the elements. The 
matrix for the models used is given below: 

Model (1) 


4 


Model (2) 

A 

are found out by talcing arithmatic average of the observed 
values. These concentrations are used to evaluate 4 matrix. 

A T A is found out and then multiplied by a scalar A t r , the 
corresponding time interval. This matrix is determined for all 
the time intervals and added to a sum matrix. The sum matrix 
is then inverted. Separately, A T matrix is evaluated and 


r 



“ C 1 C 2 0.0 0.0 
Cl - c 2 -c 2 

0.0 0.0 c 2 -c. 


Average concentrations for each time interval 
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multiplied by a vector AY^ Tdiich is nothing but the difference 
vector of two successive observed concentration values for the 
time interval, r used earlier. This vector is s umm ed up over 
all the time intervals and then .multiplied to the earlier 
inverted matrix, This directly gives the vector of parameters 
as outline in chapter 2. The algorithms are repeated for the 
model 2 also using the same programs. All the results are 
summarised in tables 6-8. 

The methods were compared for their parameter 
values, execution time required on a DBC 1090 computer, and 
degree of reliability that can be attached with the parameter 
estimates. Approximate analytical solutions were also used 
to get alternate estimates of parameters which may serve as 
references. 

A variation of linear least squares estimation 
was used for differential equations as outlined in chapter 2 
earlier , to provide the initial estimates of parameters so 
that convergence can be achieved in the program for mini- 
misation of objective function. The method worked satisfacto- 
rily for the both the models tested and gave estimates that 
were quite close to the optimal values that were found after 
the minimisation was accomplished. The method can be safoly 
used in case of exigency or lack of good computer facilities. 
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After the initial estimates were found, 
parameters were evaluated using four different algorithms 
out of which three employed minimisation and the fourth 
evaluated parameters directly from the data. It was found 
that Weighted Residual Method gives the best estimates out 
of all the four methods, tested as was reglected by the 
reliability estimates- which are nothing but the confidence 
limits for the parameters- and the residual analysis carried 
out; the parameters obtained from this method gave the least 
objective function value showing the closeness of predicted 
value of multiresponses to the corresponding observed ones. 

The method also required a negligible central processing time 
of 0.11 seconds on the DEC 1090 installation , The equations 
representing Model 1 were nonlinear in nature and did pose 
some divergence problems in the integration step and, therefore, 
it is quite beneficial to use a method that can circumvent 
this step as is done in the W eighted Residual Method . . 
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model 1 

TABLE 6a 

Initial estimate oe parameters 


Run 

Ho. 

Initial 

Cone. of 

TP A (gmol/l) 

Initial Cone* 
of EO (gmol/l) 

rate constant 
\rr for HET 
'A gmol/l min) 

rate const- 
ant for 

BHET k. 

( igmol 1 * * 4 

1 min) 

1 . 

0.476 

1.750 

0.0107 

0.0437 

2. 

1 . 590 

1.545 

0.0092 

0.0762 

6. 

1 . 600 

3.140 

0.0066 

0.1170 

7. 

1 .560 

1.910 

0.1310 

3.2240 

8. 

0.462 

1.560 

0.0072 

0. 1 670 

10. 

1.500 

1.560 

0.0102 

0.0720 


MODEL 2 
TABLE 6b 

INITIAL ESTIMATES OE PARAMETERS 


Run No. Initial cone. . Rate constant (.Jir“ ) 

of m- xylene K, K ? K, K 

( gmol/l) ' 


1. 4.000 0.0715 0.03148 0.0362 0.0962 

( 578°K) 

4.000 

(533°K) 


2 


0.0099 0.00594 0.00725 0.0160 
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MODEL 1 
TABLE 7a 


COMPARISON.. OF PARAMETER ESTIMATES OBTAINED FROM DIFFERENT 


7T 

i. 


SLGORITHMS 


Run No Measures 
of 

comparison 


Fourth. Fourth Implicit Weighted 

order order one point Residual 

R-K-G- R-K-G collocation Method 

with with with N*L*l*S* 

minimisation 


nonlinear Box 
Is minimi-complex 
zation* .minimi- 
sation 


1* 

PARAMETERS 

K 3 

0,0090 

0.0095 

0.0091 

0.0099 


K 4 

0.0781 

0.0567 

0.0783 

0.0491 

2. 

CPU TIME (Sec) 

0.67 

0.69 

0.74 

0*11 

3. 

RELIABILITY 

(k 3 ,k 4± std. 

deviation) 

0.11250 

0.08650 

0.110 

0.04376 
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TABLE 7a. 2 


MEASURES OF Fourth Order Fourth 

COMPARISON R-K-G with order 

nonlinear R-K-G 

Is minimi- with 

zation* Box 

complex 

minimi- 

sation. 


Implicit 
one point 
collocation 
with N.L.L* S 
Minimisation 


Weighted 

Residual 

Method. 



PARAMETERS 





- 

K 3 

o.ooa 

0.0076 

0.0077 

0.0082 


K 4 

0.0845 

0.0762 

0.0820 

0.0714 

Run 

N 0.2 

CPU TIME 
( SEC) 

0.67 

0.69 

0.74 

0.11 


RELIABILITY 

/ r r r r . r* n 

0.0925 

0.0872 

0.0791 

0.0670 


(K 5 ,K 4 + Std. 


Deviation) 


— • 4 

PARAMETERS 

K 3 

0.0081 

0.0077 

0.0075 

0.0074 


K 4 

0.1120 

0.0954 

0.1024 

0.1050 

Run 
No. 6 

CPU TIME 
( SEC) 

0.67 

0.69 

0.74 

0.11 


RELIABILITY 
(K 3 ,K4+ Std. 

deviation) 

0.0826 

0.0769 

0.0656 

0.. 0542 
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PARAMETERS 





K 3 

0*142 

0; 1 36 

0*132 

0.139 

e 4 

3.224 

3*046 

3.219 

3.212 

CPU TIME 

0.68 




Run ^ SE °) 

0.69 

0.75 

0.11 

N 0.7 RELIABILITY 





(K^,K^+ Std, 
deviati on) 

0.1255 

0.1134 

0.1134 

0.1044 

PARAMETERS 





K, 

3 

0.0066 

0 , 00 68 

0.0071 

0.0072 

K 4 

0.1320 

0,1350 

0. 1220 

0.1290 

~ CPU TIME 

( SBC > 

0.67 

0.69 

0.74 

0.11 

RELIABILITY 

0.0944 

0.0872 

0.0675 


(K 5 ,K 4 + Std. 
deviation) 

0.0594 




PARAMETERS 





K 3 

0.0089 

0.0086 

0.0084 

0.0082 

K 4 

0.0982 

0.0950 

0.0912 

0.0925 

Run CPU TIME 

No. 10 (SEC) 

0.67 

0.69 

0.73 

0.11 

RELIABILITY 
(K 5 ,K 4 + Std. 

0. 1122 

0.0989 

0.0745 

0.0674 


deviati on) 


MODEL 2 
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TABLE 7b 


OE PARAMETER ESTIMATES OBTAINED FRO M 
DIFFERENT ALGORIT HMS 


Run MEASURES OF 
No.1 COMPARISON 

F our xn OrderFourth 
R-rE— Gr with order 
nonlinear R-K-G 

Is minimi- with 
z at ion. Box 

Oomplex 

minimi- 

sation. 

Implicit 
one point 
collocation 
with N.L.L.S 
Minimisation 

Weighted Anal- 
Residual ytical 
Method, solu- 
tions 

1. P ARAME TERS 

( hr~ 1 ) 






K 1 

0.0761 

0.0763 

0.0726 

0.0717 

0.0721 

K 

2 

0.0362 

0.0315 

0.0293 

0.0324 

0.0344 


0.0299 

0.0321 

0.0351 

0.0371 

0.0491 


0.0594 

0.0645 

0.0772 

0.0844 

0.0665 

2 . CPU TIME 
(SEC) 

0.6600 

0 . 6800 

0.7200 

o.icoo 

0.1800 

3. RELIABILITY 

OF 

ESTIMATES. 

0.05715 

0.05421 

0.05405 

0.0484 

0.0523 

Run 

No* 2 

1 . PARAMETERS 
(Mr-1) 






K 1 

0.0091 

0.0094 

0.0092 

0.0099 

0.0097 

E 2 

0.0072 

0.0063 

0.0067 

0.0062 

0.0058 

K 3 

0.0067 

0.0077 

0.0074 

0.0074 

0.0071 

K 4 

0.0135 

0.0146 

0.0104 

0.0152 

0.0164 

2. CPU TIME 
(SEC) 

0. 6700 

0.6900 

0.7300 

0.1100 

0.1800 

RELIABLITY 

OF 

0.0324 

0.0288 

0.0307 

0.0072 

0.0084 


ESTIMATES 
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MODEL 1 


TABLE 8 


gQMgAfilSON Off_.il. S. S. W PPTT tto T , . 

BOX- COMPLEX METHOD 


IE ALGORITHM USING 


Run 

No. 

residual SUM of 

SQUARES' AS OBJECTIVE 
FUNCTION 

JE5. . QD , AS OBJECTIVE 

FUNCTION* 

1 . 

0.0754 

0.52 XI 0" 8 

2. 

0 . 0 682 

0.47 X 10" 8 

6. 

0.0722 

0.44 X 1C -8 

8. 

0.0541 

0.28 XI 0 -8 


* This objective function was used in the estimation 
of parameters while using algorithm 2. 

S 



CHAPTER 6 


GQICLH 3I0N3 AND EEC OHMSN DAT I OH S FOR FUTURE WORK 

Parameter estimation in case of ordinary 
differential equation sis studied, in detail: 

1 . The existing literature on the subject is reviewed 
and the importance of the subject is stressed. 

2. Pour different methods which cover the different ways 

in which the problem is tackled in the literature are compared 

3. To provide the kinetic investigator with a good approxi- 
mation of the parameters involved in his differential models, 
a Variation of linear least squares estimation is proposed. 

4. Weighted Residual Methods are advocated for the purpose 
for the advantage of speedy accuracy and ease of handling the 
equations. 

is 

5. Residual analysis/carried out for both the test models 
to check the adequacy of error distribution assumptions. 
Nothing concrete can be said by looking at these plots. 

6* Computer programs are developed which treat the subject 
of parameter estimation for the case of differential models 
with good accuracy and reliability. 

The field of parameter estimation and model 
discrimination has offered an enormous potential for research 
in the last couple of decades and, due to, the pioneering 
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work of some of the workers progressed very rapidly to its 
culmination. 

The future work in this direction could he 
aimed at the successful use of these techniques as a tool 
in the kinetic investigation of industrially important rea- 
ctions. In case of dynamic models, the error variance is 
not constant ( 34 ) therefore, parameter estimation and 

model discrimination based on conventional techniques may 
give erroneous results and therefore, the .algorithms : 

■fehsasfc involving integration of differential equations may be 
modified to account for changing variance of the errors. 

More so, in case of complex reaction networks where one is 
faced with a range of products and reactants and resorts 
to integral reactor for kinetic investigation due to higher 
conversions required. The experimental procedure could be 
sequentially designed to save a lot of time in experimentation 
by cutting down .the number of runs required and the tremen- 
dous human effort that is wasted in carrying . these experiments 
without any plan. Thus, a practical use of these techniques 
is more desired rather than some abstract work in this area 
which will only amount to wastage of computer time. 
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MA IN PROGRAM TO CaPhY OUT PARAMETER ESTiMrtT TOD IN MODELS 

governed by ordinary differential equations. 


PROGRAM SEGMENT !. ; INITIAL ESTIMATES 
*♦***+***********************♦**+♦***+ 

r.iis PROGRAM FI.*DS the I‘UTIAu ESTIMATES OF PAPhMETERS li 
GOV r.Ri'tfeT' rtY FxRST ORDER ORDINARY DIFFERENTIAL FOUaTtOnS 


L - M t ( i . 1 O ) , AVkC Y C i , * ) , TIME ( 1 0 J , a( i , 4 ) , XT ( 4,3 ) , XTX ( 4 , 4 ) , aT l H 
1 ) ,H.V.tV.( ),9) .SH ij x4,4 J .SUM2C4) ,Sl)MlNUiA) ,THETaU) ,wKSPCKtle J , 
2i0«»« (4 r U,AAl.;,4) #0 8(.4,4) f AVCn*C3),PATc.U5j i 

PA 1 -v CCTh;,L(I,Jj,i. = l,*t),J = l,4)/X.O,O.L,0.(j,0,0,U9 0,1.0,0,0,u»U 

X v • * # I • 0 ( u * -j # 0 a 0 f 0 # 0 f 0 9 0 f 1 m 0 / 

UmJ iltinl) ,1*1,4), J=l,4)/16*0.0/, (SUM2CI) ,1 = 1 , 4)/<**0.0/ | 
OFFutu I7=2l,i)EVICE=*o5K' f FlLfcs'FuP21.DAT # ) 

'A!' ilu j,T = 4i’ , uFV I C F 3 * DSK , FILc = ' ^OR “tO . OAT * ) 

G-.XuCTl ,* J iTiMS,lCQN, IPARA, ChTlSI : 

-'TAI.C21 ,* j U Y(I,U) ,1 = 1 , ICON) ,J = 1 , I xImE) 

REAu(21 ,*) (TlME(l) , 1=1 ,ITImE) ! 

echo sack THE VALUES read I 

-‘Hi 16X40,994) ITImE, ICON, IPARA, CATLST 

I T E ( 4 0 . 9 9 ri J f ( Y t I .0 ) . I =1 . iCiTi ) -0 = 1 .ITI MF.l 


‘IU it U*U j 77 7 J in Mt , LEUN , A1 LSI 

G-, I T £ ( 4 0 , 9 9 d 5 C CY (I ,J) , 1 = 1 , ICON ) # U = 1 , ITIME) 
JRITEtAu , 997 ) (TIMEC I ) . 1 = 1 , iTlME) 

’OR ,’iA X ( 5X , ' NO . OF MEASUR£M£nTS= * • I3//5X, # «0 


FOR mA X ( 5X , 'NO . OF M£ASUREM£nTS= * . 13 //5X , 'NO. OF SPECIES 3 U' / 
1 5 X , n 0 . OF PARhMET&RS 3 '»i 3 //sX, CATALYST CONCENTRATION FOR Xcj 
lRUu=' , Fb . 4 , grr.itOieb/litre '//) 

FORMAT xo be changed for different NO. uF species ; 

FORMAT (5X, 'THE CONCENTRATION OBSERVATION MATRIX IS '// 13 C5A,Fi A 
!))//)■ : 

FORMAT (SX, 'THE TIME IN riRS aT WHICH THE OBSERVATIONS WERE MADcl 
l/(t>X,Ft».2)) ] 

ACtuAu CALCULATIONS START FROM riERE “4 

I T I w E i = 1 T I M £ - 1 s 

CALL 5L0PE( Y ,‘TIME, RATE, ITIME, I TIME! , ICON) | 

CaLl A V R AGE ( Y , A V RGY , ICON , ITIME , I T I P£ 1) ! 

r b.i.rn,.LTr:,ufi. lot.- iina,P pn6 rsr„ -i TmF tmiv.o 


VriliU 

CALCULATIONS ARE OOnE 
Du 100 1=1# ITIME1 
Ou ID 0=1 , ICON 
AVCoN X J )=hVRGY ( J # X ) 
RATE! t J JsRATEl J,I) 
r 1 1 a, r t a i i w 


UN , 1 i T h t , l i i r* 1 1 ) 

FOR EACn TIME INTERVAL AND SUMMED UP* 


■» * * » * fcrf * I* ■ V tk v W f * / 

C 0 N IlRliE 

CALb DMATCX,AVCON,CATLSI , ICON , I PARA ) 

C A L u TP N MAT (X , XT , ICON .iPARAj 

CALL MATNUT (XT.X.XTX# I p ARA , ICON , IPARA) 

C ALL -1 AT MT C XT . PATE1 , XTY , I P ARA , ICON ) 

,-i-. -i a j i. ~ i t B * n a 


DO 20 KK=1 , IPARA 
DO 2d Jj=l, IPARA 

SUV.KAK, JJ) 3 SUMl(KK f JJ)+XTXCKE,JJ) 

S U A l ( bh ) =StM2(KK ) +XTY (KK,1 ) 

CO DTI iUE 
CONTINUE 

TYPE*, (CSUM1 Cl, J), 1=1, IPARA), 0=1, IPARA) 
PARAMETERS FOR INVERSION ROUTINE? CHANGE 
.: = IPARA ■ 

•1 = I P A R A 
I A= IPARA' 

IrsIPARA 

IC=IPARA 

IAA=IPARA 


YRA), 0=1, IPARA) „ „ 

ROUTXNE?CHAN( J E WIT ti MODEL 



IbB=IPAHA 
IFAIL=0 


«• j ^ ® '^AEf ($LM 1 , 1 A , I DEN , IB , N , M , Su* IN , 1 C , wKSPCF ' , Ah, I A A , oB , !;>£' 



tViW PARAMETER* ARK *// I. 5 X , r 1 0 .j 

G ;0 

; 1 J J T ^ 5 U 3 K 0 U T 1 N E ^ T 0 ^ F X M 0 ^ R A T E ^ 0 K ^ K F A C £ IQU 

TU^PHul £N£ SI. 0 Pfc.C 3 t, TIME, RATfc, 

' • i * c.i 1 u ■ i Y £ *• , " . ) , T ± ’!£ ( MM ) , RATE (M , V ) 

'h. IP 1 = 3,1 
■'*. 1 .1 = 1, 

.Vl\ = i:( I,J + l 3 -YCI, J) 


[ * $ * 


i =U. E(J + 1 


.> r t<(. t, J)=Ot,Ll/i-EhT 

Cun' l i HUE 
Eh , 7 OH.'. 


ibhROUT I H E 


,^ TTT , r ,„*JS,S{S 5 S^'iS 1 l J rt |' ,$ CONCENTRATIONS 

SuBKOuTlfit AVRAGE £C ,CAVRG ,fJ , M , Ml ) * ********* *********** 

REAL C(w,fO,CAVKG(N,Ml) 

>30 10 1=1 ,N 
Ou 10 J = 1 . H 1 

^AVr|GCI r j 5 = (CCI,J +1 )+CCl,J))/ 2.0 

CUHlIUUE 
RETURN 
C 1. 0 

Subroutine o m a t t v a t a . c o n , c a t . i n , i r ) 

***;T****r***:M:M************H********.* ******************* 

REAL M ATA(Tn,IK ) ,CON( 1*0 
i A 7 A ( 1 , 1 3 = ~C U in £ 1 ) 


UTACl ,2)=CON(2) 
iATAfl ,33=0,0 
i AT A (1,4 3 =0,0 
'■> ATA(2,1)=CGN(1) 

- iATAC 2 , 23 =-C 0 fn( 2 ) 
>iATAC 2 , 3 )=-COin£ 2 ) 
!ATA(2,a3=CO.N(3) 
MAT AC 3 , 13 = 0.0 
M AT A ( 3 , 23 = 0.0 
v.ATA(3,3)=C£J.%(2) 
- 1 ATAC J, 4 )=-C 0 ft( 3 ) ' 
H E T U R hi 
Hu o 


SUBROUTINE MATMTf XX f Y¥,2Zf IN. IK} 
********************************************************** 
DIMENSION XXCIM, IK) ,YK(IK) ,ZZ(IR) 

DO 10 11 = 1, IN 

zzc in=o,o 

OU 10 J U = 1 . 1 K 

ZZC II)=ZZ(II)+XXCII,JJ)*YY (JJ) 

continue 

RETURN 


£ N 0 


SUBROUTINE JlAT.MUTC AfB,C, IK, IJ.IL ) 

* * ,* **** s(t *** 4 ; * * $ * * * *** * * ** * * * * ** ******** *.* * * * * * * * * * * ****** * 



CJCJtJCiC 


r a 1 


DIMENSION A(IK,iJ),bClJ,TXi),C(lN,ll ) 

DO 1 1KK=1 ,IK 
DO I If,L=l , XL 
C(JKK,IU,)=0.U 
DO 1 LJJ=1,U 

CUr.K, lUJ=C(lF>,IIiL) tAtlKK,.TjJj *dUJJ, 1 Lu) 

D rj . V 1 .liO 

■? h I S D’.'i 


ut i.'.'u i t. • - 1', i'D>V‘rtT l A f l , iK , 10 ) 

♦ ■f************************************** 


* % 4. 4, * f 4 

i‘i' l. .j r u . , Hi) ,aT(tj,ik) 


i ’ ' * k i %) k J ** 1 jr £ t J 

V‘ X L?'ts = l,l* 

„ 1 ia^ , 1 r> J - 1 * ( if> , i J J ) 

Cw’i i 4 j Ij L 


ft !:« v i uft m 

; ; r- 


PrOuP.An St;c;,viE»*T 2* :r I ft ST ALGOhirHiM 


•Ovl!. PRuGKA i ’10 FP-D THE DPil-Ai, ESTIMATES OF PARAMETERS OF 
n SET OF OkDXEmP J£ uIFFEPt-fy i I AL E^UaTXO.mS 


— ARKATS In COi-iMON 

0 1 y E *: S 1 u U C 0 N u u ) , C OF n D C 3 ) , T I m F U 0 ) 

Cun MO ; C A T , wVAR , NP A kA , NQBS/GUeSS/CON , CZERU , TImE 


v- n a f 

CuHACirt FOUNT 

p c . al et a , ks umsq , s x epm x , x tol 

INTEGER I,IFAXL,IPRIM ,J, LI «v,bJ,LV,LW, M,MAXCAL,L, N,NF, NITER 
■LOCAL ARRAYS- 


Hh AL r\)rtC(30.4) ,FVECUO) , V(4,4) ,W(220) ,5(4) , TrtETA(4) 
r t. x^(],5 


I i'j T E 0 c,R 

C X T £R M Ab LSQFUN , LSO MO N 
DAI m C AT , r, VAR , NPAkA , NuBS , ivOoN i / 1 . 0 , 3 , 4 , 1 0 , 1 / 


* = 4 j 

: = 3t ' - I 

L.J = iO ■ | 

ij V = 4 ! 

L i E ~ 1 I 

L = 220 ! 

READ*, (CZEPU(X), 1=1 ,H VAR), ( riMEC I ) , 1=1 , NOttS) , (THETAC I) ,1=1 ,NFAR4 

1) f 0<jCl) J _ 1 3()) i 

TYPE1 I (CZERO ( i ) ”1 = 1 *, N V AR ) , ( TTME( I ) , 1 = 1 ,i*OBS) ,(TriETA(I),I = l,«PAP | 

FORMAT (IX, r INPUT: '//6X, # 1HE IMTXAl CONDITIONS C(0) aRE'// | 

i 3 ( 8 X , F 1 0 , 4 / ) / / »X , ' (HE TIrtF. VECTOR £S , //lU(«X,F7.2/)//6X, 'THt In 
UTIA b ESTIMATES OF PARAMETERS ARE V/(8X,F12,S/)//) 

FORMAT IsX'TuE' MEASURED VALUES OF THE RESPONSES ARE ' // ( 5X , 3FI P . 5 / 
1)//) 

IPRINT=1 
MAXC AL=400*N 
ETA=C . 5 
XiCb=l,0E-03 
STEPMXslO.O 
IFAIL=1 



t.| ifcj >'j >.,f 


,J > R l 

’■? ‘ f -i ‘j 

M 9 UN / 

9 v 1 '■> 

1 3. 1 1 


IFLAG=1 

CALL E04FCF (M , N , LSQFUn , LS gMGN , IPRJ.NT , MhkCAL , ETA, XTOl, STi;,R iX, 

1 THETA, FSUMSQ,FVEC,FJ AC f LJ,S, V,LV,MTe.R,NF, I J ,LIW, W,L* , IFAlL J 

IF CIFAIL.NE f 6)TYPE9999d,IFAIL 

IFCIFAIL.EQ.f )GO TO 3 

7YP£99997,FSUfnSQ 

rype.99996, ( THtTACJ) , J = 1 • NPARA ) 

T YPE9 9^95, (FVEC(I) ,1=1,30) 

Ft.. Pm A A (SX, ' V ECTOR OF RESIDUALS IS ' // ( 1 OX , FI 0 . o/ ) // ) 
uu TO l til 
TiP £9 999 1 

t\’.F.-*.<UCiX,'THh PPuGkAS IS NUT CONVERGING' J 

rufMAmx,' OoTriu: *//6x»'the final value of parameters ju>‘ 

li. *) ) 

*>•' . ^ a ( 1 >' , * 1 H £ SUM or SyUARES IS',F12.4) 

K«>? i A CC///loH ERROR EXT! TYPE, 13) 

%y J.„ fj t' 


cc 

4 


Si oROuTl'.L USEFUL UFLAG , M , .U , XC ,FV£CC , I W , Li W . W , Lw ) 

ROUTINE To c.VALuATF RESIDuAuS 

REAL tCAP(3,10).Cl(i).CDT(3),A(3,4),CZU3),CZ2(i) f CL3i.l) 

1 ,CCAb(iO) ,KU3) ,K2C3) ,X3(3) ,K4t3) ,FVeCC(30) , w (oW) ,XCtN),CZtUJ) 

it'.TRGLR I w ( l I * ) 

SCALAR ARGUMENTS 
1 uTEGcrl IFLaG.LIW , L<* , H , .* 

ARRAY ARGUMENTS 

ARRAYS IN COMMON 

REAh C0,,C30) ,CZEPG(i) ,TIMEClO) 

CO ',00.. C AT, /.VAR, NPAk A, NUBS/GUESS/COR, CZe.RO, TIME 

COMMON KOUNT 

The FOLLOWING ARE THF RUMGE-KUTIA GILL CONSTANTS 

DATA A A, 8, CC,D/U. 207 1067, 0. 29x8932, -0,707 1067, 1,707 10o7/ 

DO 4 IK-1 , AVAR 
CZR(IK)sCZEROtlN) 

00 100 IN=1 .NOBS-l 
Dl=TlMEtlN+l)-TiME(iH) 

Call do at c a ,czr ,cat,nvak, nPaRa) 

CALL DAT 4TC A , XC ,CuT , UV Art , NPARA) 

DU 5 l = i,«iVAF 
.\1U) = 0X*CDT(I) 

CZ1 UJ=CZKU) + 1 .0/2.0*Kl Cl) 

CONTINUE x 

CALL o.-lATt A,CZ1 , CAT , U VAK , nPARA) 

CALL aATMT(A,AC,CDT,N VAR, nPARA) 

DO 0 l=i,UVAR 
K2tI)=DT*CDICi) 

CZ2CI)=CZRC1)+AA*KI CI)+b*K2CI) 

CONTINUE 

CALL DMATCA.CZ2, CAT,! 5 VAR, NPARA) 

CALL MATMTC A#XC ,CDT f NVAR, NPARA) 

00 7 1=1 .NVAR 

K3U) = DT*CDtQ) „ , T 

CZ3C I)=CZRCI)+CC*K2CI)+D*K3tI) 

CONTINUE „ „ ^ 

CALL DMAT(A,CZ3,CAT,NVAft,NPARA) 

C ALL 11 ATMT C A,XC,CdT , N VAK , NPARA) 

D u d 1 = 1 , NV AH 

ClCI)=CZRCI) + l!o/6.0»CKia)+K4a)) + l,0/3.0*(B*K2Cl)+D*K3(I)) 



CCAPCI,IM)=CZRCI) 

CCARC I , IN+1 )=C t (, I ) 

CZR(I)=Cl € 1 5 

CONTINUE 

TYPS110. (CCCAPCI,J) ,1=1 »M VAR) .J=l ,N09S) 
FORMAT (IX, ’ nn-roiiT '//iy.'tkp rmarF 

13S15.S) ) 


OUTPUT V/lX, ‘THE CONCENTRATION MATRIX ISV/tbX 




V 


S U iS )-<■ . 0 
; Ju 'i a - 1 , N OSS 
Yj 0 0 = 1 , m V AP 
1 = 1 + 1 

UCC i?(U,K) .bT.Q.C )CCAPCj,lO=0.0 
Ct'Uj( I ) s4ohiCc Ar(u#tO) 

■> ^(1 j=CO-.U)-CCAb(l) 

i ,u ■; ,i' [ VO i 

IF ( r\ : 3j , SO . 1 ) Go TO U 

n, i2 1=1,30 

Sb v|S0=SU 1S0+FVFCCCI )*FVECC(I) 

TYPtiU, (FVLCCCI) ,1 = 1,30) 

T Y PEI 6 , SCMSv 

FOROATCsX, 'RESIDUALS IN FIRST ITERATION '// (5X,F10. 6/)// ) 
FuR.mAiCdX, 'SU h uF SQUARES IX FIRST ITERATION = ' ,F15.9) 

K 0 0 N T = K 0 !J N T + 1 ■ 

RETURN 

tl»0 

SUBROUTINES 

3 u S R 0 u T I N E ’ f 1 A T M 1 ( X X , Y X , X Z , I » , I K ) 

* + ******¥**#* + *******•*■*****************¥*♦****¥****+*** '-M 1 * 
DIMENSION XX(IN,IK) ,YYCIK) ,ZZ(Iu) 

|)U 1 0 11 = 1, XL 
ZiUl)si;.V 
no i \, uj=i,ik 

zzc n)= 2 ,z(ix)+xxcn,jj)*mjj) 

COL II Ulfc 

RETURN 

END 

subroutine omat(a,8,c,in,ik) 

% + ****** ****** ** 4 .****************** ******%+*********** 

0 L M E *1 S 1 0 0 A ( I w , I K ) , B C 1 N ) 

Ut,D=-R(l) 

A ( 1 , 2 ) = b ( 2 ) 

AC 1 , 3 ) =<• . U 
rCI , 4)=u.O 
4C2,l)=dtl) 

A( 2 , 2)=-B(2) 

A(2,3)=-BC2) 

A(2, U=S(3) 

A(3,l)=y.Q 

A(3 ,2)=0.G ‘ 

4(3,3)=riC2) 

AC 3 #4)=-6(3) 

RETURN 
END 

Dummy routine 






1 , W , t»rt ) . 

TYPE2, NITER, Nf 1 



y Kinm 


2 


1 1 
1 


C C 
1 1 


FOflMATUX* 'NITERS' , 15, 5X, 'NF=% 15) 

RETURIS 

E«D 

********************************** 

PROGRAM SEGMENT 3, : AiiGuRlTrtM 2 
********************************** 

PROGRAM 10 IMPLEMENT BOX-COmPLEX METHOD TO ESTIMATE PARAMETER 
I i MODELS GOVERNED BY DIFFERENTIAL EQUATIONS 
•l >r. -SIGH Y(2) ,KTL(10),SP(7) , YL(2), YH(2) , CON (4, 7) ,CZER0(4.| 
i , rr-iE(7),C(i) 

Cu- MO:, CAT,Nn,KK/GUESS/CON,CZeRO,TIME 

iK^AV DATA 

ur A m),Y(2},YL(l),YL(2),YH(l),YH(2)/u.O09,O.O4Q,Q.OO7,U.OJ 

o' l \ CkTu(I) , 1=1. iC j/S. 1,10,3.15, 10,500,0, 1 ,0/, CSPCI) ,1 = 1,7)/ 
ll. 6,1 , D 0-06,1.5,1 .5,0.02,1,0,0.0/ 

' ’ C H ij A f* i } -* j‘ r\ 

/C=i 

CAT=0.185 
Y3 4 

R c — 2 • 

Rr.Au*, CCZERU(i) ,1 = 1 ,4), (TIME (I) ,1 = 1,7) 

RE A 0 * , ( (, Cl* f« ( I , J 5 , 1 = 1 , * ) ,J=1 ,7) 

CALL CPX(LAUEl,KY,Yl, iH,KC,^T u ,4»P,Y,VAL,C,IT) 

TY PEI 1 , V AL 

FuRHAT(5X, 'THE OBJECTIVE FUNCTION = ' , 2X , E15 . 8 // ) 
riPtl,CY(I),I=l,2) 

FORMAT (5X, 'THE OPTIMUM ESTIMaTc-S OF PARAMETERS ARE '// C5X ,F1^ . 1 
1 ) 

STOP 

EuO 

SUB ROiJT I ‘IE V AGUE ( IT , V ALT ) 

********************************************************** 

TiiIS SUBROUTINE EVALUATES THE OBJECTIVE FUNCTION REQUIRED 
FUR MINIMIZATION 

********************************************************** 

R !:. A u 
ICO 1C O 

2 C 2 a ) , F v i/'UJ (Ou^vucj , v, ti-' v ^ , a i. ■* , - 
CUE iO„ CAT,HN,KK/GUESl/CON,CZERu,TIME 

PhE FOLLOWING AkE THE RUNGE-K JTI'A GILu CONSTANTS „ „ „ 

DATA AA,b,CC,L/ 0. 2 071067,0. 2 928932, **O.7G710u7,1.7u7lOo7/ 

Du 11 I IK= 1,4 

CZR(HK)=CZEPO(IIK) 

DO IlU n = l,o . 

0 T = T I M E ( M + 1 ) - T I Ft E ( N ) 

CALL DMAT(A.CZR,CAT,HN,KK) 

CALL M ATM! ( A, YT ,CilT , Nn ,RK) 

DO 2 1=1,4 , 

Kl ( 1)=DT*CDT( I) 

CZt(I)=CZR(I)+1.0/2.0*KI(I) 

CiJN 1’IiJLF 

CALL UMATCA.CZl ,CAT,NN,KK) 

CALL MATMTCA,YT,CDT,Nn,KK) 

DU 4 1=1,4 

K.2(I)=DT*CDT(I) „ 

CZ2(I)=CZft(I)+AA*Kl(I) +B*K2 (I ) 

CONTINUE 

CALL DMAT(A,CZ2,CAT,NN,KK) 


***************************************************** 

YT(2) ,COn( 4,7) ,CZERO (4) ,TIrtE(7),CCAP(4,7),Cl(4) ,CSTAR( 4 1 * 
D,A(4,2}.e£i(4) ,CZ2(4) ,CZ3C4),K1(*),K2U) .K3(4) ,K4(4) ,LCAu 
, FVECC ( 2b ) , SUM SO (26 ) ,CZR ( 4 ) , S ( 4 , 4 ) , R5D ( 4 , 7 ) 



t d 4 


B 

1 wO 


IJLll 
1 1 1 


:cc 


10 


c 

c 


CALL MATMT(A,YT,CDT,N«,iCK3 
DU 6 i = l,4 
K3CI )=DT*CDT(1) 

CZJU J=CZRlI)+CC*ic2U j+0*K3U) 

C u m X X 0 0 E 

CALL) Ui'-iH T t A , C tj 3 , C AT , N rt ,kK ) 

C-M.u . ; A '1 lT(A,YT,Ci»T,P..»,KK) 

‘A ■ A 1=1,4 

•<XX Us^TfCDTU) 

U J I =C ARC J ) + i . u/ 0 . 0* i K 1 £ 1 3 +K4 U ) 3 + 1 
it'd, )-C.Pll) 

~ o c 1 , +1 )=C1 ( X) 

A t I J =C 1(1) 


V- i 

Cc 

Ov 

Cu 

n. 

ih.i 

n 


G/j.0*(«*K2Ci)+n*K3(X ) ) 


a. .Or, 
111 i : 
tj ill l\ : 

S u s v , i . 


1,4 

1,4 


mi j= 1,7 


, u 


V 


•CCARtN , J) 

, 3 ) 


E6D X K , J J = COu C i\ 

I i’- A . = h S J ( X , J 3 * r. O T C * 

Su.” = $u ! i + Tc,hm 

S C I , K 3 = S U m 

crmt i.MUti 

■i « 0 E n= ■* 

CALL DETER (3 , UETS , N ORDER ) 

Vrtl,T--ao.0**l3j*OFTS 

RETURN 

END 

sunKoaTirjts 

So A K 0 U T I . j £ M A IK X ( A X , Y Y , Z Z , I w , 1 K ) 

dimension XX (IN, IK) ,YX(XK) ,ZZ(I.O 
DU 10 11 = 1, IN 

zzan=o,o 

00 10 J J = 1 , 1 K 

ZZl II) = ZZCIX) + XA( II,JJ)*YY(JJ3 

CQWTIUUE 
R £ T U R i-i 
e.xD 

SUBROUTINE OMAT(A f B f C,lN,IK) 

t * ****** *****¥***)|:*4.***:f*+$*$t*f********#*¥**4'***** + ** 

DIMENSION A(IN f IK),bQN3 
TER MSB C2)*C/C2,0*btl)tB(3)) 

TfiRMl=2.0*BU)*TEFM 
T£Rrt2=6(3)*TERM 
A ( 1 , 1 3 =~TERMl 
AC 1 ,2 3 = 0.0 
A C 2 , 1 ) = -TERMl 
A 1 2 , 2 3 = -T£RM2 
AC3,13=TERM1 
A ( 3 ,2 3 = -TERM2 
AC4, 13=0.0 
AC4 f 23=TERM2 

RETURN . 

End 

II i S M fw Y ' 'R n 1 1 T T M F % 

iii^*^****^* ****** ^* ************************ ******** 

SUBROUTINE TRACECLAb, XTER3 

RETURN 

END 



2 

U> 

n 

*v« 

r* 

\4 

u 

2u 

21 


W 

22 


30 

31 

32 


• SuBROOTIN E « F 1 R C ( L A <3, IT £ R ) 
RETURN 
Erin 

SUB? GuTi IMVAP(Yi) 

m •*£ sign y x c 2 ) 


Ot.TUR.* 

) 

Su 3 k n j c i 
■ * i *■ u * r o v 
RET OR < 


: c f static . , jc.‘ ) 
C .V C 0 C ! * 1 


Sv( ^ lu r t •:« iiE 1 P.K C A , jT£R>t , N ) 


* '4 + * f * * 


I'O FI. 


TnF 


p,j _ 

i \ ;i 1 i.i ’ A 1 ' r '■ ) 

.1 -A. !, iT B’i t I v r ' t R L C“>.vD£F3ATI0N 

0 l 0 ~’ = 1 

*3 = - *■ i 

HO iLI'VXi-JATlO^ i-,-1 TImES 
0 L> B 0 I = 1 t N l 
11=1 + 1 

Fl\u PIVOTAL fcLCHtM 
B IG=0 ,(• 

I R G >v = 1 

lCuu=l 
)U 10 J = 1 t N 
00 1C K=1,N 
Ab=ABS(A(JfK)) 

IF ( d I u~ AB ) 2 » 10 , 1 0 

BlG=Aii 

IKOw=J 

1C0L=X 

CUN'l Irilit 

U UROa- 1121,21, 11 
EXCHANGE RD« XPuW WITH 1, 

MULTIPLY b Y -1 

D'iER Y=~oT£RFi ' 

DO 20 ■)=[,•'! 

Tt *'P=A(1 ,0) 

(UI, J 3 =A C IF.a^ 

UIRQ*,G)=TeMP ^ 

1 f ( IC jL-I) 31 » 3 l f 22 t , , T 

CuLuMN EXCHANGE ICOL aI'IH I 

MULTIPLY DTERs bY -1 

Ol£EM=-DTERM 

iiu 3t J-Iii# 

rhNp=A(j,i) ^ 

M,J»I)“A(G»XCuL) 

MO, ICODsTEMP _ 

DTEKM=DTERM* A ( T , 1 3 
IFCAII, I))32,V6 f 32 
RATlOsl./AClr U 
i)U 40 K=II,N 

DIVIDE ROW ELEMENTS PY PIV01 
AU,K)=A(I A K)*RATIp ^ 

ELIFiluATS CuLuMri ELeMeTTS 
DO '-tO J=I I » >4 . 

TERO=A( J » 1 5 *A ( I , K ) 

ACJ,K) = ft(0,K)-TtR<4 


uRuER m-MFIX 



3 3 

4C 

5 0 


k D- 

'L, 

07 


0IF=l0000,0*AaS(A( J # K))-AdS(TeP;4) 

IFCDIFJ33, 33,40 

A(J,K)=0.0 

CONTINUE 
C G r I N U £ 

pT£KM=DIERM*A ( N , N ) 
r Y D t‘- 60 , uTERM 

••'..K .M’CiX , 'DETERMINANT IS'/5X,E15„B) 

r r u r .-I 
TVPc.^7 

E v j x i A i ( 1 -j a f 1 3 n n ET E R M I N A N TsO ) 

Re. t-UH.ii 


H> 


20 

30 


40 

ItH) 


101 


w 


? o a U A ^ S | G St. ^ T ^ 4 . * S AG 0 U R I T rj M 3 
;.aT , PROGRAM TO ImPUEME.-<T WEIGHTED RESIDUAL 4i£TrfOo FOR 

* vt *♦***♦ U***i* f ** i** ****** i ******* ******^* ATA0 * S 

Au C (4 , 7 ) ,CAVkG (4,6) .CO, ( 4 ) , DELCOM C 4) , MATA (4 * 2 ) • MAT AT (2 » 4 ) 

1 / ':* AT AT A ( 2 , 2 ) , ATOELC C 2 ) , SUM i ( i , 2 ) , SUM2 ( 2 ) * THETA ( 2 ) , TIRE (7 ) , 0t.i.T 0 
2b) ,0tLC(4,t>) , Ah(<*, 4 ) , VRSPC G(28),6R(4,2),i r )c.N(2,2), SUM1 IN (2 , 2 j 
9 AT A CAT/0.l85/.Nf,,*M.lA,IB,TC,IAA £ I3d/2,2,2,2,2,4,4/, (CIDE.HA,4 
U,X = l,2),J = l,2m.0,0.0,0.0,l.u/,IFAiL/0/ * I 

R t A D * r M t M f K 

READ*. ( (C(I,J) ,X = l,iO,J = l ,M), (TIMECl) ,1 = 1 ,M) 

•i 1 Sri~ 1 

CALL AVRAGE (C ,C AVkG , W , M , M 1 ) 

C A Lb DELTA ( C , X I HE , DeLC , DEuT , N , M , M 1 ) 

00 10u 11=1 , Ml 

DO 10 Jj=l , N 

CON t JU ) =C AV RG ( J J , 1 1 ) 

CALL DESIGN (CUM ,MATA, C4T,N,N) 

CALL im MAT ( MATA, t'lAXAX,*,(0 

CALL RA1 MUT(MATAT,M*TA r MATAl’A,K ,N,K) 

ox=acuTan 

00 20 JU — 1 ,6 
Du 20 Kk=l , K 

MATAT.4C J J , KA)=MATATA( JJ , KK) *DT 
0 0 3 0 J u = 1 , N 
OELCOUCJJ)=0£LCCJJ,II) 

CALL mATMTCMATAT , DELCQN , ATDELC ,K, >») 

00 40 Kk=l * K 

SuMl UK . KK) = SuM 1 C kK , Kk ) tMATAT A( i\K , KK) 

SUM2 UK ) =SUM2(Kk) + ATDELC (KK) 

CONTINUE 

CONTINUE 

CALL F04AEFCSUM1, IA,IDE«,IB r NN,MM,SUMllN,IC,WKSPCE,AA, IAA, 

1*1 FAIL) 

TYPE** CCSUMIINU.J), I»| ,2),J=1,2) 

CALL INVRTCSUwl ,SUM1IN,K) 

CALL MATMT ( SUM 1 IN, SUM2, THETA, K,K) 

UPc.101 , (THETACI) , 1=1, K) , „ ,, 

FORMAT C5X, 'THE VECTOR OF PARAMETERS: '//(5X,F3. 5/) ) 

STOP 

********************************************************* 

SUBROUTINES START FROM HERE ,.. e 

SUBROUTINE TO CALCULATE THE AVERAGE CONCENTRATIONS 
SUBROUTINE AVRAGECC,CAVRG,N,M,M1) 

\ 



f i : 


10 


1 u 


REAL C(N,S'i) ,CAVR<UM,M1) 

DO 10 1=1, N 
DO 10 J=1 , Ml 

CAVRGtT, J)=CCCI,J+1 )+C(X,J>)/2.0 

Cu*i n.NUE 
R E UiRJl 
S.iD 

4 * * ; 

Su-^0UT1'1^ fO CALCULATE TnE DIFFERENCE SECTORS 
• . *KQtJ?l .c. DELTA (C , FlME , DELC ,L£ lT , N , M , M l ) 


' * 4 * U »-/ * . U i M 1 W, ? l JL nr* fOLUL # U C iJ X m » « # - *1 , 

C( * ) > DtLC ( , H 1 j , TIME ( V ) , DE lT ( M i ) 

■ ■ j t < * I = 1 , " i 

DO !■ . 1 = 1 , x 

0 .ICC i,u) = (CCl»J+i)- r, CI,J)) 

E i.tCvO = :i,'f.( kJ + i)-Ti‘ c(J) 

L U ' < J 1 1 1 U E 
Rul uR* 

Lri.i) 

* * 4 


so 5 n.?uTiLE 'in evaluate me uEsion matrix 'V 

H'-i LOuTj, ,*«, 0LSIb\(Cuf!,MATA,CAi,N,X5 
r» L Alt CO IN ) ,MATAC b , K) 

’• r.R •/.sCO t.l2)*CAT/C2,0*CO.\(l)tCON(3)) 

X 1 =2 . 0*COL ( 1 ) *TERM 
X2=C0'iO)+TlRm 
•IATaC 1 ,1 }=-Xl 
■ AT A C 1 , 2 ) =0 . 0 
• AT A C 2 » 1 ) =-A l 
1 AT A ( 2 , 2 ) =-X2 
!AtA(i,n=xi 
*> AT AC 3 , 2 ) =-X2 

•inTAC't,l)=O.U 

MaThU,2)=X2 

RET URL 
EuO 

¥ 4 *¥* 1 ^$¥**** 1 c ¥¥¥*¥¥¥*¥*t¥**¥*¥*****¥*¥¥¥¥**-¥***¥*¥***¥* 

subrouting m atmt cxx , yy , zz , i n * ik ) 

D 1 M g 0 STUN XXQN, IK), YYCIK), 2211(0 
DU 10 11=1, 

22(111=0,0 
DO 10 J J = 1 , 1 K 

ZZ C 1 1 ) = ZZ C 1 1 ) +XX (I I , J J ) *YY C J J ) 

CONTINUE 

RETURN 

EfJO 

♦X****************************************************** 

SURkOUTIUE MAIMUTCA,R,C,IK, 1 J,Il) 

DIMENSION A(I<s,IJ) ,B(IJ,IL) ,C(IK,IL) 

DU i IKK=l,IK 
Du 1 ILL=1 , XL 
C C IKK, ILL J =0 . 0 
00 1 IJJ=1,1J 

C(IKK,ILL)=C(IKK,ItiL)+ACIKK,IJJ)*i3(IJJ,ILL) 

CONI’ I LUg 

return 

******** ********************** ******** 

SUBROUTINE TRNMAT ( A , AT , IK , I J) 

DIMENSION ACIK,IJ),ATCIJ,IK) 

00 1 IJU=1,IJ 



;t jr jtif j Ci cjcj 


DO 1 IKK=1,IK 
RETuRu 

C* i* i) 

INVER |S .A DIAGONAL MATRIX*************** 

“ 1 K ^TC.-.Mr, IMTTri^) 1 


.;•• ■ 10 1 = 1 , K 
i ,' 1 - M 0 = 1 »>- 
i f C i m i'd \f m Jf J 
;. °a : (i f ,i) = 

*•<!•> i'-.i It- 

.. i r ci, jj = 

wu ■ 1 1 i i.t E 




5 ?S . <iSTI«A TES 


plpiilitlsn 

svHLUAfeD K FRo!i r A D | E f H SF i !io§ux^4 S n iM i^r9i^ 1I ^^^^^ RAMEii: ' RiJ 

Rt.hu KSIGMB J 1 3 ' 4 J ' H 1 ' * F C 4 ' 4 5 * A C 4 , 4 ) , PM AT (4,45, DPGPM (4.3J 

82^j*^oSsy^!!{?j,^ n " AnM ^ irjR TnE ru * 

|;:?AY c o^A V “ k ' wai,4/l -' J '-* / 

OATA 1 */ f r-iflr • <t "^MATRIX TURNS OUT TO BE IDENTITY MATHjLX 

^ 'iSiii ??§52 ‘ il£*?>stix „< 4/ 

j, , 4 " ^ f f £ (cov,um.j) f J— i f ? j * x * 1 1 - j ) 

P*?g| AC5X T,lE Er ' P0R VARIATE MATRIX IS'//(5X, JF15.5/)//) 

FOR MATC5X, 'OUTPUT: ' J 
00 10 1=1,3 
gMATU,I) = 1.0 
CZR(I)=CZERO(i) 

CONTINUE 

CALu TR.MMAT(GMAT,GTRA^S,M,LJ 

ou ioo a = 1,9 

■Hi 25 1=1,3 
00 25 0=1,3 
00 25 K=1 ,4 
AflATd, J)=0.0 
Bo AT ( I, K) = 0,0 
no AT ( I, K) = 0.0 
CONTINUE 

n X = T I M E ( N + 1 ) - T X M E ( N ) 

CALL RUNGE (CZR, THETA, DT,C1) 


CALL* JACXCC1, THETA, AMAT) 

Call jacacci , theta, bmatj ; 

CALL RKGMCAmAT ,BVAT,DT,DMATJ 
CALL TRNMAT(DMAT,DTRA«NS,M,K) 

CmLL tATMiJTCDTRANS, GTRAnS,SPMAT,K,L,M) 

CALu aATmuTCC.MAT,DMAT,SMAT,m,L,K) j 

CALL UA TAUT (.SP.MAT , COVMA I « DrLPM , X , M . L) * 

CALL aAT^UT(DPGPM,SmAT,H1mAT,K,L,K) 

Du X' 1 = 1 ,K 
>u 20 J=1,K 

A '.AT ( l , J ) = H f, A T ( i , j 3 +Ht MAT ( I , J) 

: j . 4 T I t U E 

■V J«, T = l,3 . ; 

CL'HI)=C1CU 

CuLTIhIJE ; 

CALL VR T( -i v < VT,Ps'-'»M‘,K,K) 

Ou 5=0 1 = 1 , i j 

SOH = SUM + P>iAlCl,i) 

A i fL -1 A b = S 0 R T ( S U M J • i 

T.iS D'EGkEES Of' FREEDOM ARE = o. THEREFORE t alstrimition value 
TAKE A is t = 2 . 5 7 1 { 

KSI„Mb=2.57l*3IGMAB 

TYPES 1 I 

Ou 09 1=1 , NPAR A J 

T YPE50 , THETAC 1 3 , KSIGMb j 

FORMAT CSX, Fi2. A, 7X,F12.L/) 

COW TINGS 

FuROAT CSX, 'RELIABILITY uF TrE PARAMETERS : ' //9X ,' PARAMETER 8* j 
l AKul ABILITY'//) 

STOP ! 

SfiO . . ■ ■ 1 

SUBROUTINE RUnG£(CZr7xH£TA~uT,C1 ) ” | 


t. 7071067/ 


REAL Kl(3) ,K2C3) ,K3(3) ,K4(3) ,CZRC3) ,C1(3) ,XC3,4) ,CDTC 3 3 ,TriE TA U ) 
1 ,CZ1 (3) ,CZ2C3) .CZ3C3) ! 

CUM MOW /COMST/CAT,NVAR,i*PARA 
DATA A, B,C,D/0. 207 1067, 0.29289 32, -0,707 1067 
CALu D M A T C X , C ZR , C AT , M V A R , N P A R A 3 
CALu M ATMT CX, THETA, COT, 0 VAR, WPARA) 

no 2 1=1,3 
K1CI)=OT*COTCI) 

CZUI)=CZRU)+0.5*KlCl) 

CONTINUE 

CALL DM AT (X ,CZ1 ,CAT,N VAR , nPARA) 

CALu MATMTCX, THETA, CDT, AVAR, WPARA) 

00 4 1=1,3 

K2CI)=0T*CDTCI) v _ , . 

CZ2U)=CZRU) + A*KlQ)tB*K2U) 

CONTINUE 

CALL 0MAT(X,CZ2,CAT,UVAR,NPARA) 

CALL AATMTCX, THETA, COT, NVAR,^PARA) 

DO 6 1=1,3 , 

K3 C 1 3 =DT*COT Cl) _ % 

CZ3(I)=CZR(1)+C*K2(I)4D*K3(I) 

CONTINUE .-.o^ 

CALL DMATCX,CZ3 / CAT,NVAR,NPARA) 

CALL HATMTCa, THETA, CD T,NVAR,NPARA) 

DO 8 1=1,3 
K4 ( 1 3=DT*CDT C 1 3 



CH I)=CZR(I) + 1.0/6.Q*(Kl(n+K4(I)m. 0/3.0* (8*K2U)+0*K3Q)J 
Continue 

RETURN 

S')HkOuTINc,S~--— 

3 u A K Out I ME ATM X ( XX , YY , ZZ f I n , IK ) 

,H ,o!C J XX ( i N , IK ) f Y Y (IK ) , 2Z (In ) 

, 1 \J 11 = 1# IN 
o Z C i l ) =0 , 0 

5u 10 = 

'■* A C 1I) = Z7.1U)+Xji(II,JJ)*YY(J.J;> 

Cu-U’InOE 

R c, I u R ;'j 
i «i) 

3v 3R0u|I ‘it u* £ AT l A , R ,r , T x , 1 K ) 

v I ■? *-: ;■* j 

( i ,l) = -.i(l) 

AU ,2} = uC2) 

A Q , 3)=o.U 
A a, 4) = 0.0 
A(2,u=«icn 
A ( 2 , 2 }=-8 (2 ) 

A C 2 , 3 ) = ~8 C 2 ) 

A(2,4)=6C3) 

A ( 3 , 1 ) =0 . 0 
413,25=0.0 
AC3,3)=b(2) 

4(3,4)=~S(3) 

RETURN 


EmO 


subroutine jacxici , theta, amatj 


REAL Cl ( 3 ) , THETA C 4) , AH A 1(3,3) 
COMMON /CONST/CAT, N VAR,. NPARA 
ABATCi,l)=-THETA(n 
AM AT Cl ,2)=TrlETA(2) 

Art AT ( i , 3)=0 , u 
Ak A X ( 2 , 1 ) =T H E X A ( 1 ) 

A MAT (2,2)=- ( THETA ( 2 ) + THETA ( 3 ) ) 
Am AT ( 2 , 33 =TtiEX‘A (4) 

A;iA X ( i , 1 ) =0 . U 
A*nA T ( 3 , 2) = TmE f A ( 3) 

Am A T ( 3 , 3 ) =-ThbTA ( 4) 

return 

E N 0 

SUBROUTINE J ACA (C 1 , TriF.TA , 8MAT ) 


REAi, Cl ( 3 ) , T8r,TA( 4) , RmAT ( 3 , 4 ) 
COMMON /CONST/CAT, NVAR,nPARA 
BM AT (1,1)=-C 1(1) 

BMATCl , 2)=C1 (2) 

BmATC 1 , 3)=G .0 
UMAX (1,4) =0.0 
BmAT(2 , 1 )=C1 ( 1 ) 

BMATC2 , 2 ) =-Cl (2) 

BMAX (2 , 3 ) =-Cl (2) 
8MAT(2,4)=C1C3) 

UMAX (3,1) =0.0 





8 MATC 3 , 23 = 0.0 
BHATC3,33=Cl (23 
8 MATC 3 , 4 )=~Cn 3 ) 

RETURN 

EuD 


SUBROUTINE’ RKGM(AMAT,8MAT,DT,OMAT3 


"‘EAL AM AT ( 3 f 3 ) , dM AT t 3 , 4 ) , OH AT ( 3 , 4 ) . K1 C 3 , U , <\ 2 C 3 , 4 ) , K 3 l i , i j 
14 ) ,P'uiCC 3,4 3 ,FuAC 1(3,43, FU*C 2 ( 3 , 4 3 ,FUNC 3(3,4 J,ADmAi ( 3,4 3 
Du 1 0 1 = 1,3 

■') U i 0 ij — 1 | 4 


<)u 2 *j * 
A U r H 1’ 1 I. 



, vl )-A r \. Ai ( I , J 3 + AM AT Ci , K ) A I ( K , J 1 

J3 = At • ATU,J3+ttMAT(I,J) 

= 0 I * r U. C 1 1 , J 3 


K** t 3 < 


EO--CI ( I , U) = r'U. CU f J)H,5*Kl(I,J) 

k 2 C i. » u 3 = u 1 * I U hi C l ( i , U 3 I 

FuuC2(.[,.3j=HK-Ca,J3+A*NHI,J3+o*K2CI,J3 

K3U,J)sOX*FUaC2(I # J) ! 

Ku'.C3CI. , JJ=rTi..ClIf J)+C**2(T,JJ+D*K3( I, -J3 

K 4 (l,J) = OT*PU;,C 3 (i.U) j 

DuAT( l, J3=FUNCCi,J)+1.0/6.0*(M (I, J)+i\4(T, J) 3 + 1 . O/ 3 ,0* Cb*K2 ( i 
l+D*K3(l,J3) 1 

CONTINUE 


RETURN 


End 


SuBKOuT I N E >s AT HUT ( A , B , C , I K , I J , I u 3 


DIMENSION A ( IK , ID) , b( I J , IE) ,C ( I»\ , xL3 
DO 1 IKK=1,1K 
DO l ILL=1 , IL 
C I IKK , ILL 3=0 .0 
00 1 IJJ=l,Id 

C C I KK , ILL) =C ( IKK , ILl) t A ( IKK , I J J } *ii ( i<J J » ILL ) 

continue 

RETURN 

End 


SUBROUTINE 1 Phf* AT t A , AT , IK , I J 3 


DIMENSION AUK,IJ) ,AT(tJ,lK) 

00 I IJJslflJ 
Du 1 IKK=1,IK 
ATIIOJ, IKK3=A(IKK,IJJ3 

Continue 

RETURN 

SUBROUTINE IN VRT(A, AINVRS,M«,KK) 
REAL A(MM,KK3 , AINVRSCmM , KK) 

no 10 i=i , mm 

DO 10 J=1 # MM 
IFU.mE-JJGO TO 20 
AINVRSCI, J)=l ,0/ACl, J3 

GO TO 10 
AINVRSC 1,03=0.0 




CONTINUE 

RETURN 

end 

*********** * *********** *** **************************** 

PROGRAM TO PLOT RESIDUALS OLTaUEO AFTER ESTIMATION 
* * * * ************************************************** 

D i " n. JGin.'f T IMc ( 1 0 ) , ¥ ( 1 0 ) , lSuRT (10) , COn (4,10) 

• b A 0 * f '4 M 0 U E Li 
•■'U 10 1 = 1 , . iODFL 
READ*, HORS,,, VAR 

0;.AO*, ( (I iF(J) ,J=l,fiObS) , t (COi(J,K),J=l , NVAA),Ksl # NOBS) 

DO 3 o J=1 , ,VAk ■ 

Du 20 K=1 , MuiS 

V )=CCb.(0,r ) 

Cii;> 1 lOUt*. 

DSTEP\=60 
• Li T E ? 1 = t» D 
IrAiL=0 

CALL G01AGF (Y,TIME,NCi>S,ISOKT, JSTEPX, *S1EPY, If AIL) ' 

IT P t * , I F A I L 
CONTINUE 
Cu./iMnUt 
.STOP 
f N D 

********************** ********** 

PR 0 G R A M S E G M E NT 5 . ^ : A L G ORITrtH^J^ 

PROGRAM TO FIND THE OPTIMAL ESTIMATES 0* °ARA»ETEPd OF 
A SET OF ORDINARY DIFFERENTIAL EQUATIONS <JSUG 
ORE POINT COLLOCATION METHOD 

************************************************************* 

— ARRAYS IN COMMON 

DIMENSION C(28) ,CZEkO(4J ,TImE(7) 

COMMON C AT # N N , K A. / GUESS./ C , C ZE RO r T I ME 
REAL ETA , FSUMSO , STEPMX , XTOL 

INTEGER I , IFAIL, IPBIN f ,d ,LI n ,LJ ,L V , Ld , M, MAXC AL,L,»» ,RF , NITER 
--LOCAL ARRAYS 

REAL FJAC(28,2),FVEC(28),V(2,2),vH125),SC2),X(2) 

Integer i*u) 

EXTERNAL LSgFuN.LSOMOX 
DATA C AT ,i<N,KK/0. 105,4,2/ 

4 = 2 
M = 2 8 
LU = 28 


L V = 2 j 

L 1H=1 ■ j 

i,., = 125 , j 

READ*,(CZERGU),I = l,4},CTIMEU),I = l,7),(XU),i = l,2) \ 

READ*, (CU), 1 = 1, 28) „ | 

TYPE1, (CZERO(I),I = l,4) | (TIM£(i),I = l r 7) { (XU) ,I = i,2) , 

FORMAT (IX, # INPUT: V/6X, THE INlTiAL CONDITIONS C(0) APc. // 

18X,(F10.4/)//6X,*ThE TIME &ECTOR I S V/8X , ( F7 . 2/ ) / /oX , 'THE I^lTi} 
1 AL ESTIMATE OF PARAMETERS IS'//8X, (F12.5/)//} j 

FORMA^ (6X' ThE = MEASURED VALUES OF THE RESPONSES ARE V/BX , ( 4F1& .5/| 
1 )//) 

IPRINT=1 
MAXCAL=400*N 
ETA =0.5 



r c! 


i 

t;i 

y* - ! 4 


J l 


9 j j i 
u j i iM 

U)i - 


1 1 L 


9M 


XT0L=1 .QE-03 
3T£PMX=10.0 
IFAIL=1 
IFLAG=l 

CALL EO 4FCF(H,N.LSQF Un ,LS\)MGN , IPRIN I , MAXCrtL, E l A , X TOl, STr.PoX , * . 
1 FSU 130, F VEC , F J AC , LJ , 3. V , liV .Ul TER , NF , I W , L I W , X , L* , IFAT u ) 

IF (IFAIL.NE.<mYP£oMtf,lFAIu 
If UFAIu.EQ.l )GU TO T 
a I 9097, F 30 m SO 
fYPc.^y t i.ib,(ACu),J=l,2j 
go to lilt 


program is not CONVERGING') 

GuTPUi : ' // 6 X , ' TrtF FINAL /ALUS Of PARAAFXFrS li'i 


'£99941 

FuP FP-'lC 1 X, ' i H c. 

FuP . Ai ClX, ' 
i / a 6 a , f a . 4 ) j 

Fun FA Tax, 'TOE 3 a. OF r^OAP^S I3*,F12.4) 

FORi'iAT(///lon FkFOP. EXI l TYPE, 13) 

3 1: D p ■ ■ 

■ar 

t a 

S L 3 P. 0 U T 1 N E L S w F u N C I F L A G , 'I , X , X C , V V E C C . l W , L I '<* , M , L a ) 

*****4-*************t*f* ********************************** 

routine to evaluate pestouals 031 uG semi-implicit punge-kotta 

'’.ETuOu 

Ai-'PAY A f\ G u M EM S ~ - 

r< l*. A L FVECC C 28 ) . t. ( uW) , XC C N } 

I. .1 EGER IwliuT.a 
SCALAR APGU.MF.jT3 — 

1 i . [ C.G3I' I F L A G , L 1 W , L , rt , N 
ARRAYS 1.4 COMMON-** 

Rt-.Au C ( 28 ) ,CZeFG ( A) , 'j ME ( 7 ) 

LOCAL ARRAYS 

REAL CZRC4) ,Fu! (4) , AM AT c 4, 4) , i ATR IX ( 4 , 4 ) , IDEM ( 4 f 4) , Kl (4) , *1 A T A L V 
14,4) ,F(4) ,CCnPl4,7) ,CCAL(2R) , NKSPCEQO) ,AA(5,5) , RB ( 5 , 5 } ,C1 (a 
COMMON C A T , ia , K e / G n E S 3 / C , C Z E R 0 a I M £ 

D.VCA UlOElUl ,0) ,J = i ,4) ,1 = 1 ,4)/l6*0.0/ 

• ** 1 llC 1 1=1, .J 

czru)=czergq) 

I »)£N C I , I ) = 1 .0 
00 100 I E = 1,6 

I)TsTImE(IN + 1)-T1ME(IK) 

C All u F U NC (XC/C ZR , F U«,u N , KfQ 
CALL jACXCCZR,XC,AMAT) 

00 ‘>0 1 = 1, N .■» 

Du 90 J=1,NL 
HATRIXC I,J)=0,0 

MATRIX! I, J) = lLEN(I, J)-0.5*DI*AMATCI f J) 

Continue 

N L L = 4 
i M M = 4 

1 A= 4 
Io=4 

I C = 4 j 

I AA=5 
108=5 
IFAIL=1 

CALL F04 AEF (MATRIX, I A, IDEN, 18, NNN,MMM, MATIN V, 1C, WKSPC£,AA, IAMBS 

1,168,1 FAIL) 

00 80 1 = 1 ,NN 



80 


U' 


F C I )”0T*FUN (I) 

CALL MATMT(MATINV,F,K1 ,NN,fLM) 

DO 70 1=1, HI* 

Ci (I)=CZR(I)+Al tl) 

CCARU, lfO=CZ«(l) 

CCAF(l,lX + l)=CHI) 

C2RC1)=C1 Cl) 

Cu iX I on jr 
C«: m T I M 3 £ 

K s v , 

Ov t> vi 0 = 1 ,7 

i'i 60 1 = 1,0:. 

f\ = !".+ ! 

0(;u,u)=ccahU,j) 

FUCCC*) =C C i\ ) -CC A L C K ) 

c u ■. r i o u i, 

Rc,i‘LRx 
F:, Q 

S»iOR 'luT ). 'it e. S - — - - — — — --------- 

Su3K0uTI.'hi kA iMT (XX,Yi ,ZZ,Iu,iK) 

Hi UUolu t X U m,lN) ,YU i<) , 2 , 1 ( 10 ) 

00 10 11 = 1 ,1'' 

7/2(1 X)=0,o 
00 10 JJsl.AK 

22 C LI)=ZZCXi) + XXCn,Jj)*YX (JO) 

CQ-4xiLue 

RETURN 

i.<0 

•>0 -1 iY RiuU’lNF 

SuatiOUT 10c, uSLMOh ( M , N , F V ECC , Fj ACC , LJC , S . IGRaOc, . NITER, «F 

******** ***** + *****♦***** is********;*** ***********♦*#**♦***♦¥♦** 

X , N , L X ) 

?YPE2 f NITER, Of 

FORMA hiX, 'NlTErt=' , 15,5X, *NF=* ,15) I 

R l.U T ii R 

SUBr,nuTI v uC FUriCtThETArCUNCffnM^.K) 

^f ^****^********#***** 

« ? i^u C0f4CCM)/fn^(.W) # TrfEIACK) 

CO i **10 N CATfAiN|KK % 

X 1 = 2.0* jO-iLlACl) *CuNC(l)*CO«C(2)*CAT/C2.0*CONCCi)+CONCl3)) 

X2- IHETAC 2 ) *CuflC C 2) *CGHC ( i) *CAT/ (2 .0*COuC C 1 ) +CONC ( 3) ) 

FUN ( 1 ) =-Xl 

FU)U2)=-X1-X2 

FU,U3)=Xl-X2 

FUN(4)=X2 

RETURN 

EoO 

S tj ft ROUTINE JACXCC1, t’HETA.AMAT) 
******************************** 

REAL Cl (4) , THETA ( 2) , AM AT (4,4) 

XI = CH2)Sci?35/t2.*Cia)tCH3))**2 

53s8i»}«??S«lii*gia cu3))-2 

X4=CH3)/C2,*CH1)+CH3)) 

A.MATCl,n = C-2.*THETACl)*CAT*XU 
AM AT C 1 , 2 ) = C -2 . *TBETA(1 5 *C AT+X2 ) 

AOAT ( 1 , 3 ) ~ C2 , *THE'f AC 1 ) *C AT*X3 ) 





AHATC 1 , 4 )=Q ,0 ! 

AMATC2,n = C-2.0*THErAU)*CAr*Xl+2.*i , Hc,TAC2)*CnT*Xl) 

AHAi'(2,2) = (”2,*T’Ht,TA(l)*CAT*X2~TH£TA(2)*C\r*X4) 

A;UTC2, J)=C2 t *THETA(l)*CAT*A3~2.0*TttETA(2}*CAT*x3 ) 

AtiA 1 C 2 r 4 5 =0 . 0 i 

An A r C 3 , 1 ) = { 2 . ♦THF*1* A C 1) *C AT*Xl +2 . * i HETA C 2 ) *C AT*X 1 ) 

A.iAVU,2) = (2,*TnEmn*CAT*\2-TrtEm2) *C AI *a4) ; 

.*..!AV(3,3 5 = (-2.*THETACn«‘CAT*Xj-2.*THrrAl2J*CA'l , »x3j : 

\ A 1(3,4) =0 . 0 

A Ai C *,l) = (-2.*i'HETA(2)*CAT + Xl) 

A . . A r ( '1 , 2 ) = ( THETA ( 2 ) *C AT ♦ <4 ) 

A •••;,*( 4 r .ps2.C-»TnETA(2)»CAr*X3 

A t i A 1 ( '% , 4 ) — U • 0 

RtvluR.. | 

S i i * A 0 u T L ’> £. 1 VR'i'CA.Al f, V R S , M A K ) 

H & A « . ACM •" , K K ) , A I h V R 6 ( xM ,K.<) 

00 1 0 1=1 f WM ' f 

J u 10 J = i , K ,-l . I 

lr ( L.«K..))G.j TO 20 . 1 

\ i ’» y R3 ( I , J ) = 1 .0/A(I,J) l 

go in io I 

A I ' i v R 5 ( i , J ) = 0 , 0 i 

CO ‘-IT I ;hje - 

RETURN . .. | 

EXO ’ . ■ i 
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